Document [original]

Evaluation of Highly Viscous Liquid

Mixing Based on Particle Tracking

Dissertation

Submitted to the

Department of Chemistry

Faculty of Natural Sciences

University of Paderborn, Germany

for the degree

Doctor of Natural Science (Dr. rer. nat.)

Bhaskar Bandarapu

Adviser:

Prof. Dr.-Ing. H.-J. Warnecke

Institute of Technical Chemistry and

Chemical Process Engineering

Department of Chemistry

University of Paderborn

Co-adviser:

Prof. Dr. rer. nat. D. Bothe

Chair for Mathematics

Center for Computational Engineering Science (CCES)

RWTH Aachen

Date of submission: 10.09.2008

Date of examination: 17.10.2008

To my

Family, Teachers and Friends

Acknowledgements

It was a great opportunity and valuable experience for me to carry out this research

work in the group of Technical Chemistry and Chemical Process Engineering at the

University of Paderborn, Germany. I would like to thank many people who directly

or indirectly, technically or morally supported me during the course of this thesis.

First of all, I would like to sincerely express my deep gratitude to my adviser Prof.

Dr.-Ing. H.-J. Warnecke, who accepted me as his Ph.D. student and provided this

challenging topic. His constant support and guidance all along my stay in Paderborn

helped me to complete this work successfully.

I am greatly indebted to Prof. Dr. D. Bothe, Chair for Mathematics, CCES,

RWTH Aachen, who agreed as a co-adviser of this dissertation and for his invaluable

suggestions and critical comments. I am very fortunate to work with him and to attend

his interesting lectures to improve my mathematical skills. Without his constant

guidance and encouragement, I would never have been able to complete this thesis. I

hope this work would partially fulfils his expectations.

I would like to thank Prof. Dr. K. Huber and Prof. i.R. Dr. H. C. Marsmann

for accepting as the members of my Ph.D. examination committee. Financial sup-

port from BASF AG Ludwigshafen and German Academic Exchange Service (DAAD)

together with STIBET scholarships, is greatly acknowledged. Here, I would like to

mention and thank Dr. Bobert for looking after the administration work in this regard.

I particulary thank Martin Schmidtke and Houman Shirzadi for friendly atmo-

sphere and constant technical/non-technical discussions at our ’Multi-Kulti-Büro’. I

thank Dr. Raach and Mr. Maletic of Thermal Process Engineering group, for allowing

me to use their computational facilities whenever needed. I would like to acknowledge

Tobias Maxisch for his help with the mapping method. Many thanks goes to our sys-

tem administrator Mr.Voigt for his timely help whenever I encounter problems related

to computer and softwares. I would like to address and thank all the members of the

group, former and present, for the pleasant and supporting atmosphere.

I personally thank my friend Dr. V. K. Badam for regular discussions on technical

and non-technical issues. I would like to thank my teachers and friends who encour-

aged and supported me in many difficult times. Finally, my deepest gratitude to my

parents, sisters and brothers who always stood behind my decisions and encouraged

me constantly throughout my studies. I would like to appreciate their patience and

support, which allowed me to stay abroad for a long period of time.

September 2008 Bhaskar Bandarapu

List of Figures

2.1 Illustration of continuum hypothesis, where the density ρis calculated

from molecular mass ∆mwithin a given volume ∆V. . . . . . . . . . . 8

2.2 Description of time derivatives. Assume that a fluid element is moving

along a curve Cthrough points 1 and 2. . . . . . . . . . . . . . . . . . 11

2.3 Schematic representation of the Reynolds transport theorem. . . . . . . 12

3.1 Illustration of control volume approach with three dimensional dis-

cretization notation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2 Schematic of the simple reference geometry considered for rotating and

sliding mesh analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.3 Section plots of velocity vectors after one rotation of the impeller in

moving mesh (left) and rotating reference frame (right) simulations. . . 34

3.4 Contour section plots of velocity magnitudes at different time intervals

in both moving (left) and rotating reference frame (right) simulations. . 35

3.5 Computational sensor points to extract quantitative information for

better comparison of the velocity magnitudes. . . . . . . . . . . . . . . 36

3.6 Quantitative comparison of velocity magnitudes computed both in mov-

ing mesh and rotating reference frame simulations. . . . . . . . . . . . . 37

3.7 Schematic of (a) the Couette-flow geometry and (b) the computational

grid which contains 207,000 hexahedral control volumes. . . . . . . . . 38

3.8 Comparison of theoretically calculated torque with torque computed

using Fluent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

vii

3.9 Comparison of theoretically calculated torque with torque computed

using Star-CD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.10 Comparison of theoretically calculated torque with torque computed

using Star-CD for highly viscous liquid. . . . . . . . . . . . . . . . . . . 41

4.1 Schematic representation of anchor mixer, where all the dimensions are

in mm (Geometrical configuration is received from BASF AG, Lud-

wigshafen). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.2 Cell occupation of liquid phase (0 < αl< 1) in a surrounding gas phase

(αl=0) based on the VOF-method. . . . . . . . . . . . . . . . . . . . . 48

4.3 Computational grid which consists of approximately 300,000 tetrahedral

control volumes generated using ICEM-CFD tool. Mesh with magenta

colour (below) correspond to liquid phase and the one with aqua colour

(above) correspond to gas phase. . . . . . . . . . . . . . . . . . . . . . 49

4.4 Computational grid together with the mixer geometry at the interface

between liquid and gas at an instant of time. . . . . . . . . . . . . . . . 54

4.5 Computational grid at the interface between liquid and gas at an instant

of time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.6 Representation of concentration of scalar variable in a mixer vessel

stirred with anchor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.7 Contour plots of the velocity magnitudes in a mixed vessel stirred with

anchor. High velocity magnitudes can be seen arround the anchor arms. 56

4.8 Pressure distribution in anchor mixer. High pressure results at the front

of the stirrer and low pressure behind the stirrer. . . . . . . . . . . . . 57

4.9 Time evolution of torque in anchor mixer. . . . . . . . . . . . . . . . . 58

4.10 Time evolution of torque in anchor mixer for various viscosities. . . . . 59

4.11 Influence of solver settings on numerical diffusion. 1st order discretiza-

tion (left) and 2nd order discretization (right). . . . . . . . . . . . . . . 60

4.12 Influence of computational mesh refinement on numerical diffusion.

Coarse mesh (left) and fine mesh (right). . . . . . . . . . . . . . . . . . 60

4.13 Contour plot of number concentrations (right) obtained from the num-

ber of tracked particles inside compartments. . . . . . . . . . . . . . . . 61

4.14 Decay of normalised variance under application of the Baker map. Ef-

fect of numbers of particles (P) and compartments (C): 100P, 64C (red);

800P, 4C (blue); 3200P, 64C (black). . . . . . . . . . . . . . . . . . . . 62

4.15 Illustration of Baker map . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.16 Typical particles trajectories in Lid driven cavity. . . . . . . . . . . . . 63

4.17 Schematic representation of anchor mixer with hybrid (combination of

tetra and hexa) mesh. The mesh contains about 100000 control volumes

and is generated using Star-Design. . . . . . . . . . . . . . . . . . . . . 64

4.18 Lagrangian particles trajectories in anchor mixer. . . . . . . . . . . . . 64

4.19 Schematic representation of anchor mixer with dimensional configura-

tions (in mm) and it is divided into 4 compartments vertically for par-

ticle tracking. A blob of particles also can be seen in the figure. . . . . 65

4.20 Schematic representation of anchor mixer with division of 16 compart-

ments in radial and circumferential cross section. . . . . . . . . . . . . 65

4.21 Time evolution of normalised variance for anchor mixer. Effect of num-

bers of particles (P) while compartments (C) kept constant: 100P, 64C

(red); 1600P, 64C (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.22 Time evolution of normalised variance for anchor mixer. Effect of num-

bers of particles (P) while compartments (C) kept constant: 100P, 64C

(red); 1600P, 64C (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.23 Time evolution of normalised variance for anchor mixer. 1600 Particles

are initially placed in one compartment (blue), multiple compartments

(red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.24 Time evolution of normalised variance for anchor mixer. Effect of num-

bers of compartments (C) while particles (P) kept constant: 1600P, 4C

(red); 1600P, 64C (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.25 Time evolution of normalised variance for anchor mixer with 1600 par-

ticles and 64 compartments. . . . . . . . . . . . . . . . . . . . . . . . . 74

5.1 Single shaft kokneader with oscillating and rotating mechanism; 1. Bar-

rel, 2. Kneading pin, 3. Kneading element, 4. Shaft, 5. Oscillation and

6. Rotation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.2 Schematic of the self-cleaning mechanism of kneader by movements of

the pins [83]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5.3 Two shaft intermeshing Reacom 60L kneader reactor. . . . . . . . . . . 79

5.4 Schematic of (a) unit cell of a specific kneader and (b) the top view of

kneader element (CAD data from BASF AG, Ludwigshafen). . . . . . . 81

5.5 Horizontal section plots of velocity fields. (a) Contour plot and (b)

vector plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.6 Time evolution of normalised variance for kneader element with 1600

particles and 64 compartments. . . . . . . . . . . . . . . . . . . . . . . 83

5.7 Schematic of the mapping matrix method: stretching of fluid element. . 86

5.8 Schematic of the mapping matrix method: particle distribution. . . . . 87

5.9 Time evolution of variance in a kneader element: Comparison of simu-

lated data with data obtained from mapping matrix. . . . . . . . . . . 89

5.10 Time evolution of normalised variance in a kneader element. Compar-

ison of simulated data with data obtained from mapping matrix after

every rotation of the impeller. . . . . . . . . . . . . . . . . . . . . . . . 90

5.11 Time evolution of normalised variance in a kneader element. Compar-

ison of simulated data with data obtained from mapping matrix after

every one-third rotation of the impeller. . . . . . . . . . . . . . . . . . . 90

5.12 Time evolution of normalised variance in a kneader element. Effect of

mapping matrices [52]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

5.13 Time evolution of normalised variance in a kneader element. Compari-

son of simulated data with data obtained from mapping matrix [52]. . . 92

List of Tables

3.1 Comparison of torques calculated theoretically and using Fluent. . . . . 39

3.2 Comparison of torques calculated theoretically and using Star-CD. . . . 40

3.3 Comparison of torques calculated theoretically and using Star-CD for

highly viscous liquid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

List of Symbols

Abbreviations

CDS central differencing scheme

CFD computational fluid dynamics

CICSAM comprehensive interface capturing scheme for arbitrary meshes

CS control surface

CV control volume

FD finite difference

FE finite element

FV finite volume

QUICK quadratic upwind interpolation for convective kinematics

RTD residence time distribution

UDS upwind differencing scheme

VOF volume of fluid

Roman Letters

aacceleration, m/s2

Aarea, m2

bsource term in discrete equation

cspecie concentration

CuCourant number

ddiameter, m

xiii

Ddiffusion coefficient, m2/s

Dij rate of deformation in tensor notation

Drate of deformation tensor

Eenergy, J

fface

fbody force, N

Fforce, N

gacceleration due to gravity, m/s2

g1, g2, g3, g4nodal co-ordinates

Hhigher order terms

Iintensity of segregation

Jdiffusion flux, mol/m2s

kface number

Lcharacteristic length of macroscopic flow, m

mmolecular mass, kg

nnumber of particles

nunit vector normal to the surface

NeNewton number

ppressure, N/m2

PePeclet number

Ppower, W

r1, r2radii, m

rposition vector, m

Ssource term

ttime, s

Tij stress in tensor notation

Trate of deformation stress tensor

uvelocity, m/s

uggrid velocity, m/s

Vvolume, m3

x, y, z Cartesian co-ordinates

Greek Symbols

αVOF coefficient

δKronecker symbol

given fraction or threshold

ηkinematic viscosity

γconstant

Γconvective term in Navier-Stokes equation

λmolecular mean free path, m

µviscosity, kg/ms

Ωangular velocity, rpm

φintensive property

Φextensive property

ρdensity, kg/m3

σvariance

τtorque, N-m

Subscripts

ax axial

Agit agitation

Diss dissipation

Kin kinematic

max maximum

mix mixing

0initial

xvi

Contents

1 Introduction and Aim of the Work 1

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Motivation and Contribution . . . . . . . . . . . . . . . . . . . . . . . . 3

1.3 Organisation of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 Mathematical Formulation 7

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Governing Equations of Fluid Flow . . . . . . . . . . . . . . . . . . . . 9

2.2.1 Basic Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.2 Continuity Equation . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.3 Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2.4 Species Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3 Numerical Methodology 19

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.2.1 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2.2 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . 25

3.2.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 26

3.3 Rotating and Moving Meshes . . . . . . . . . . . . . . . . . . . . . . . 27

3.3.1 Rotating Reference Frame . . . . . . . . . . . . . . . . . . . . . 28

3.3.2 Mesh Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

xvii

xviii CONTENTS

3.3.3 Simple Reference Geometry . . . . . . . . . . . . . . . . . . . . 33

3.4 Torque Computation in a Couette-Flow . . . . . . . . . . . . . . . . . . 37

3.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.4.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.4.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . 42

4 CFD-Analysis of Anchor Mixer 45

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.2 Analysis of Partly Filled Anchor Mixer . . . . . . . . . . . . . . . . . . 46

4.2.1 Volume of Fluid (VOF) Method . . . . . . . . . . . . . . . . . . 47

4.2.2 Problem Specification and Numerical Approach . . . . . . . . . 48

4.2.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . 53

4.3 Particle Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.4 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.4.1 Intensity of Segregation . . . . . . . . . . . . . . . . . . . . . . 68

4.4.2 Mixing Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5 Mixing in a Kneader Element 75

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

5.2.1 Flow Computation . . . . . . . . . . . . . . . . . . . . . . . . . 79

5.2.2 Particle Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . 80

5.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.3.1 Flow Patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.3.2 Mixing Quality . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.4 Mapping Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5.4.2 Basic formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 85

5.4.3 Application to a Kneader Element . . . . . . . . . . . . . . . . . 88

6 Conclusions and Outlook 95

Chapter 1

Introduction and Aim of the Work

1.1 Introduction

Industrial fluid mixing in stirred vessels is one of the basic and important requirements

of most production systems in bio-chemical, food, pharmaceuticals, polymer, pulp

and paper, waste water treatment and petroleum etc industries. The term ‘mixing‘

refers to an operation or process which tend to reduce inhomogenities or gradients

in composition, properties or temperature of the material bulk. Fluid mixing is the

motion and/or contacting of a single or multiphase process with a liquid continuous

phase to achieve a desired process result. Practically, every plant will contain some

sort of tank where mixing is carried out by the relative movement of material between

the various parts of the whole mass using suitable equipment for the specific operation

or process. There are several ways to provide mixing action in a vessel, but this thesis

is primarily concerned with the impeller type mixers in stirred vessels.

A study of the mixing process includes several basic considerations. The first is

the effect of the vessel on the mixing process. Vessel geometry, dimensions, and the

structure may dictate mixer selection and mixing performance. The next is the mixing

impeller(s) to be used for a given specific process or operation and it would be a real

problem if every new application requires a new impeller design. Thus impellers are

usually made in an homologous and in geometrically similar series. Further, the flow

2CHAPTER 1. INTRODUCTION AND AIM OF THE WORK

patterns generated by impellers can be divided into two basic types as axial- and

radial-flow impellers. Axial-flow impellers produce flow parallel to the impeller shaft

where as the radial-flow impellers discharge fluid to the vessel wall in a horizontal or

radial direction.

Although mixing is one of the most ancient technological practices, there are still

notable achievements to be done. The importance of future mixing and various steps

to follow for better understanding of it is clearly demonstrated in a review article [36].

Since mixing is primarily due to the relative motion within the material to be mixed,

the efficient design of the mixer gain much appreciation. Most of the available articles

based the results on optimising the designs of the mixers focused on the flow field,

power consumption and mixing time etc. Ample experimental and numerical data for

the design and optimisation of mixers can be found in [55] and additional references

given there. However, there is no systematic method for selecting and designing the

mixers. Further, the majority of these investigations are concentrated mainly in the

turbulent flow regime although many industrial mixing operations are carried out in

laminar flow condition.

In laminar mixing, inertial terms in the equations of motion for incompressible

fluids are either not present or not important because of their magnitude. The vis-

cous terms in the equations of motion dominate the flow behaviour since viscosity is

very large in such fluids. Additionally, rapid mixing or homogeneity cannot be easily

obtained for highly viscous liquids since there are no turbulent eddies which enhance

mixing by increasing the rotational speed of the agitator. Apart from large viscosities,

the low diffusion coefficients are also responsible for laminar mixing to be an area-

controlled process. The increase in interfacial area during mixing provides additional

area, over which diffusion can act so that mixing can be accomplished. Therefore, the

impellers used in laminar mixing are usually full-tank, close-clearance impellers. This

means the diameter of such impellers approaches the diameter of the vessel to bring

motion to the entire vessel volume. The present thesis contributes for the mixing of

highly viscous liquids in close-clearance stirred vessels under laminar flow condition.

1.2. MOTIVATION AND CONTRIBUTION 3

The preceeding section introduces the motivation as well as contribution of the current

studies.

1.2 Motivation and Contribution

In the past, extensive experiments were carried out by several researchers to realise the

effect of various parameters on mixing efficiency (quality) for different mixing devices

at variable process conditions and material properties. However, for small changes

in the design parameters of the mixers or for the different materials to be studied in

the mixer, it is too laborious and expensive to perform experiments each time. On

the other hand, technological progress in computers from the past two decades made

it easier to explore and apply in engineering sciences. In particular, computer based

simulation method "‘computational fluid dynamics (CFD)"’ has become a widely used

tool for analysing, optimising and supporting the design of mixing processes [6,8,70].

In general, several numerical methods exist in the literature to assess the quality of

mixing and some of the widely known among them are numerical tracer experiments,

Lagrangian particle tracking simulations and entropy based calculations. Numerical

tracer experiments are extensively used methods to quantify mixing and plenty of

information is reported in the literature [55]. Although this method is of increasing

importance as a means to analyse mixing processes, there is a principle problem which

becomes a severe obstacle in case of highly viscous liquid mixing or, more generally in

case of high tracer schmidt numbers. In this situation, a sufficiently accurate solution

of the species equation can be spoiled by the effect of so-called numerical diffusion [21].

A way to avoid this problem is to replace the continuous tracer concentration by a

number concentration obtained from Lagrangian (i.e. inertia free) particles. This

approach does not suffer from artificial diffusion, since the position of tracer particles

can be resolved with sub-grid-scale accuracy and the velocity field at these particle

positions can be obtained by interpolation from its values at grid points. Please refer

to [32,35,50,61,64] for further reading on Lagrangian particle based mixing calculations

4CHAPTER 1. INTRODUCTION AND AIM OF THE WORK

for characterisation of different mixers. In addition, entropy based characterisation of

mixing also gained much attention in recent years [61, 82]. The definition of entropy

is given as a measure of disorder or degree of distributive mixing.

In the present work, theoretical method of calculating intensity of segregation,

intensity of mixing and mixing time based on Lagrangian particle tracking is pre-

sented. The approach is as follows. The total computational domain is divided into

smaller compartments (sub-volumes) and inertia free particles are placed initially in

one compartment, say. During the process of particle tracking the resulting number

concentrations are recorded and allow for computation of evolution of its variance.

Based on the statistical measures, we provide an answer to the number of particles

and compartments needed for a reliable assessment of the mixing quality. Further,

this method is evaluated using the numerical investigation of mixing in a vessel stirred

with an anchor type impeller as well as a specific kneader element operated under

laminar flow condition for highly viscous liquids. Mixing times were calculated for

both mixers based on evolution of intensity of segregation. Finally, a mapping matrix

method is elaborated to evaluate the quality of mixing. This mapping method em-

ploys a transition matrix, which describes how many particles are advected from one

compartment to the other compartment in a particular period of time. With the aid

of this transition matrix one can compute variance evolutions and mixing times using

vector multiplications with significantly less computational effort.

1.3 Organisation of the Thesis

The contents of the each chapter are briefly described in this section. Where relevant,

the key problems addressed in the chapters are highlighted, and regarded them as

the important contributions of the work. The present thesis contains six chapters,

including an introduction (Chapter 1) and conclusions (Chapter 6). Since this work

contributes for the modelling and simulation of mixing in stirred vessels, it first starts

with the mathematical formulation. Chapter 2 presents an overview on the basic

1.3. ORGANISATION OF THE THESIS 5

equations of fluid flow. These governing equations are developed based on continuum

hypothesis, according to which the physical quantities are considered to be the con-

tinuous functions of space and time. The resultant system of equations are presented

both in integral form as well as in differential form.

In Chapter 3 detailed information on numerical techniques was provided. A finite

volume based commercial CFD - tool is used in the current studies and hence a com-

prehensive description is supplied on this method. A spatial and temporal discretiza-

tion schemes are elaborated to solve the governing transport equations. Simulation

of mixing in stirred vessels is always a challenging task as it involves complicated de-

sign mechanism. Different available numerical methods (e.g. rotating reference frame,

moving mesh) to model flow in mixing vessels are emphasized and employed on simple

reference geometry. Finally, torque computation in a simple Couette-flow between two

rotating cylinders is discussed and compared with the analytical results.

The numerical analysis in one of the ancient yet widely used mixing apparatus,

anchor mixer is addressed in Chapter 4. A review on the relevant literature is also

given there. Detailed flow field computations are carried out and torque is computed

for partly filled vessel stirred with anchor type impeller. The torque computation

procedure is given and the power number is explained. With ever increasing demand

for the estimation of mixing quality based on the numerical methods, there exists a

severe obstacle in the form of numerical diffusion when numerical tracer experiments

are carried out especially for highly viscous liquids. A method of avoiding this problem

is addressed in this chapter based on Lagrangian particle tracking. Additionally, an

answer for the minimum number of particles as well as number of compartments is

provided with the aid of statistical measures and finally evaluation of mixing in anchor

mixer is performed using particle tracking.

Chapter 5 extends the Lagrangian particle tracking simulation to investigate mix-

ing in a specific kneader element. A summary on the relevant literature for kneader,

single and double screw extruder is given there. Further, a mapping matrix method is

presented which employs a transition matrix to describe the advection of particle in a

6CHAPTER 1. INTRODUCTION AND AIM OF THE WORK

particular period of time. This transition matrix allows to estimate the mixing quality

with significantly less computational efforts. The comparison of simulation data with

mapping matrix is presented. A summary and outlook of the thesis is presented in

Chapter 6.

Chapter 2

Mathematical Formulation

2.1 Introduction

The first step in dealing with the problems related to hydro or fluid dynamics, is to

understand the physical background and describe them using mathematical equations.

These equations are derived by postulating the continuum i.e., continuous hypothesis

which allows to use the representation of physical quantities such as velocity, pressure,

density and temperature. These quantities are considered to be continuous functions

of space in the fluid and of time. Of course, a similar hypothesis is made in the

mechanics of solids, and the two subjects together are often designated as continuum

mechanics. If λis the molecular mean free path and Lis the characteristic length of

macroscopic flow, then the ratio λ/Lis called Knudsen number. The assumption of

continuum hypothesis is valid at macroscopic level (if λ/L1) to represent transport

phenomena in gases, liqiuds and solids. However, contiunuum hypothesis breaks down

on the molecular or microscopic level (if λ/L≤1) where each fluid is descrete with its

properties fluctuating randomly and are not in local equilibrium. This is due to the

fact that the molecules in such a process do not get sufficient time and space to come

to local equilibrium.

Nevertheless, the vast majority of phenomena encountered in fluid mechanics fall

well within the continuum domain. This situation is illustrated in Fig. 2.1, where the

8CHAPTER 2. MATHEMATICAL FORMULATION

Microscopic

Uncertainity

Macroscopic

Uncertainity

DVVD

Microscopic

Uncertainity

Macroscopic

Uncertainity

DVVD

Figure 2.1: Illustration of continuum hypothesis, where the density ρis calculated

from molecular mass ∆mwithin a given volume ∆V.

density ρis calculated from molecular mass ∆mwithin a given volume ∆V. If ∆V

is very large, ρis affected by the inhomogeneities in the fluid itself. As ∆Vbecomes

smaller, and almost uniform density is reached independent of ∆V. In continuum

approximation, the point density is defined as that value of ρwhich occurs at the

smallest magnitude of ∆V, before statistical fluctuations becomes significant. This

limiting volume is ∆V∗, below which molecular variations may be important and

above which aggregate variations may be important. Thus, the density ρof a fluid

now becomes:

ρ= lim

∆V→∆V∗

∆m

∆V.(2.1)

In the present work, the dimensions of interest are very large compared to microscopic

level, and hence it is assumed that the behaviour of the fluid in terms of macroscopic

properties. In this case, the continuum hypothesis is valid as a statistical average

of corresponding properties of a large numbers of molecules. Now, on this basis it is

possible to establish equations governing the motion of the fluid which are independant

as far as their form is concerned.

2.2. GOVERNING EQUATIONS OF FLUID FLOW 9

2.2 Governing Equations of Fluid Flow

In this section, the mathematical statements of two fundamental physical principles

are developed upon which most of the fluid mechanics is based and they are conser-

vation of mass and momentum. Additionally, conservation of species equation is also

formulated. The derivations of governing transport equations of fluid flow given here

are based on the lecture notes of [13] and other references [4, 21, 81]. Before delving

into the elaborations of the equations, it is worth a while to mention some of the

advantages and disadvantages of mathematical modeling in academic or industrial re-

search. Inspite of the high initial investment in mathematical modeling, especially as

relates to experiments:

• serves as an alternative to experiments

• is cheaper and faster

• can be used in planning experiments

• is used in parameter investigation to answer "‘what if"’ questions

• provides information about certain transport phenomena that are otherwise dif-

ficult to measure experimentally.

2.2.1 Basic Concepts

All the fluid mechanics is governed by conservation of mass, momentum, energy and

other constitutive equations. These can be mathematically expressed in many ways

and are briefly discussed here.

Differential Element Approach

Consider an infinitesimally small fluid element in the flow field with a differential

volume ∆V. The fluid element may be fixed in space and fluid is moving through it

or it may be moving along the stream line with a velocity equal to the flow velocity

10 CHAPTER 2. MATHEMATICAL FORMULATION

at each point. The fundamental physical principles applied to this small fluid element

leads directly in terms of deformation rate in partial differential equation form. The

partial differential equations obtained directly from the fluid element fixed in space

are in conservation form whereas the partial differential equations obtained directly

from the moving fluid element are in nonconservation form.

Control Volume Approach

Consider an arbitrary closed volume drawn within a finite region of the general flow

field. This volume defines a control volume CV , and a closed surface which bounds

the control volume is the control surface CS. The control volume may be fixed in

space with the fluid moving through it or the control volume may be moving with the

fluid such that the same fluid particles are always inside it. The fundamental physical

principles applied to the fluid inside the control volume and the fluid crossing the

control surface gives the transport equations in integral form. Of course, on the ma-

nipulation of these integral forms of transport equations generates partial differential

equations. The integral form of the equations obtained from the finite control volume

fixed in space are in conservation form whereas, the equations obtained directly from

the finite control volume moving with the fluid are in nonconservation form. The

control volume approach is considered to be mathematically more rigorous as it does

not assume the solution to be continuous before hand. However, both the methods of

formulation gives the final form of the equations which are independant of the method

of derivation. In this chapter we will focus on integral form (i.e., control volume ap-

proach) of the equations since most of the commercial CFD codes are based on finite

volume method.

Time Derivatives

Let us consider a fluid element moving along a curve Cas shown in Fig. 2.2. When

this fluid element moves from point 1 to 2, one of its velocity components changes

from u1at t1to u2at t2. The time rate of change of velocity of the fluid element as t2

2.2. GOVERNING EQUATIONS OF FLUID FLOW 11

Figure 2.2: Description of time derivatives. Assume that a fluid element is moving

along a curve Cthrough points 1 and 2.

approaches t1becomes:

lim

t2→t1

u2−u1

t2−t1

=du

dt .(2.2)

Here, du

dt is known as the total derivative and is defined as the time rate of change

of velocity of the given fluid element as it moves from point 1 to 2. This kind of

description of fluid flow is called Lagrangian flow analysis. It is convenient to use fixed

coordinate system and measure the velocity with respect to the coordinate system as

it is difficult to follow the fluid element. Therefore, the difference in velocities between

successive fluid elements which arrive at a point 1, seperated by time difference ∆t, is

given by:

lim

∆t→0

u1(t+ ∆t)−u1(t)

∆t=∂u

∂t .(2.3)

Here, ∂u

∂t is known as the partial derivative and is defined as the time rate of change

of velocity at a fixed point 1. This kind of description of fluid flow is called Eulerian

flow analysis. Thus, du

dt and ∂u

∂t are physically and numerically different quantities.

Similarly:

lim

∆x→0

u1(x+ ∆x)−u1(t)

∆x=∂u

∂x ,(2.4)

where ∂u

∂x is the partial derivative and is defined as the time rate of change of velocity

at a given instant of time t. Furthermore, both Lagrangian and Eulerian flow analysis

can be related as follows. Since u=u(x, t):

dt =∂u

∂t +∂u

∂x

dt +∂u

∂y

dt +∂u

∂z

dt ,(2.5)

12 CHAPTER 2. MATHEMATICAL FORMULATION

ControlVolume ControlSurface

Rateofinput

Rateofoutput

RateofGeneration

andDissipation

CV CS

ControlVolume ControlSurface

Rateofinput

Rateofoutput

RateofGeneration

andDissipation

CV CS

Figure 2.3: Schematic representation of the Reynolds transport theorem.

where (x, y, z) = xis the coordinate of the fluid element being followed at instant t.

Therefore:

dt =∂u

∂t +u∂u

∂x +v∂u

∂y +w∂u

∂z ,(2.6)

can be written and (u, v, w) = uis specified. Additionally, the vector operator ∇in

three dimensional cartesian coordiante system is defined as:

∇=i∂

∂x +j∂

∂y +k∂

∂z .(2.7)

Finally, we write:

dt =∂u

∂t +u· ∇u . (2.8)

In the above equation du

dt is the total derivative, which is the time rate of change of

velocity of a moving fluid element; ∂u

∂t is the local derivative, which is the time rate of

change of velocity at a fixed point; (u·∇u) is the convective derivative also commonly

known as substantial derivative, which is the time rate of change of velocity due to

the movement of fluid element from one location to another location in the flow field.

Reynolds Transport Theorem

Consider a system with a control volume CV fixed in space and enclosed with a control

surface CS as shown in Fig. 2.3. Reynolds transport theorem states that the rate of

change of an extensive property Φwithin the system is equal to sum of the rate of

2.2. GOVERNING EQUATIONS OF FLUID FLOW 13

change of Φwithin the control volume and the net rate of Φthrough the control

surface. Let φbe an intensive property and is defined as Φper unit mass of the fluid.

Since Φis defined with respect to a fixed mass of the fluid, the conservation laws (mass

and momentum) are also applicable to φ. Then, the relation between Φand φcan be

written as:

Φ = Z

ρ φ dV (2.9)

where ρis the density of the fluid. Let ube the velocity and nbe a unit normal drawn

outward on CS. The Reynolds transport theorem can be stated as:

Rate of change

of Φin time

=Net outflow of Φat

the control surface

Rate of change of Φ

by generation and dis-

sipation

Mathematically this means:

dΦ

dt =Z

ρ φ u·ndS +Z

∂

∂t (ρ φ)dV . (2.10)

The above equation is often called the Reynolds transport equation or control volume

equation.

2.2.2 Continuity Equation

Continuity equation is based on the law of conservation of mass, according to which the

mass of a closed system remains constant regardless of the processes acting inside the

system. This means that matter changes from one form to another form but neither

created nor destroyed. For example, for any chemical reaction in a closed system, the

mass of the reactants must be equal to the mass of the products.

The integral form of the continuity equation is obtained directly from Reynold’s

transport theorem [Eq. 2.10], by setting dΦ

dt =0 and φ=1:

∂ρ

∂t dV +Z

ρu·ndS = 0 .(2.11)

14 CHAPTER 2. MATHEMATICAL FORMULATION

In vector calculus, Gauss divergence thoerem (which states that the flux of a vector

field on a surface is equal to the volume integral of the divergence on the region inside

the surface) is used to transform the surface integral into volume integral. Therefore,

using Gauss divergence theorem, surface integral in Eq. 2.11 can be expressed as a

volume integral:

ρu·ndS =Z

∇ · ρudV. (2.12)

Replacing surface integral in Eq. 2.11 using Eq. 2.12, gives:

(∂ρ

∂t +∇ · ρu)dV = 0 ,(2.13)

for any arbitrary control volume CV . Allowing the control volume to become in-

finitesinally small, leads to a differential coordinate free form of the continuity equa-

tion:

∂ρ

∂t +∇ · ρu= 0 .(2.14)

In the above, the velocity vector urepresents three velocity components in cartesian

coordinate system. Expressions in cylindrical and spherical coordinate systems for

the continuity equation can be found in many text books [11,57]. Continuity equation

represented in tensor form:

∂ρ

∂t +∂(ρ ui)

∂xi

= 0 ,(2.15)

where xi(i=1,2,3) are the Cartesian coordinates and uiare the corresponding Carte-

sian components of the velocity vector u. In many applications involving liquids, the

density ρis constant and the flow is said to be incompressible. For such flows, Eq. 2.14

can be simplified as:

∇ · u= 0 .(2.16)

The flows representing the above equation are called divergence free.

2.2. GOVERNING EQUATIONS OF FLUID FLOW 15

2.2.3 Momentum Equation

The momentum of a point mass mwith velocity uis given simply by mu. The

momentum equation is based on conservation of momentum, according to which the

rate of change of (linear) momentum of a fluid element is equal to the sum of the

forces acting on it, which is nothing but Newton’s second law of motion.

Rate of change of momentum = Sum of all relevant forces.

Therefore:

dt (mu) = For alternatively, F=mafor m= const ,

where ais an acceleration. The integral form of the momentum equation is obtained

directly from Reynold’s transport theorem [Eq. 2.10], by substituting dΦ

dt =PFand

φ=u:

∂

∂t (ρu)dV +Z

ρu u ·ndS =XF.(2.17)

Here, Frepresents all the relevant forces and they are:

•surface forces: which act directly on the surface of the fluid element (pressure,

shear and normal stresses, surface tension etc.);

•body forces: which act directly on the volumetric mass of the fluid element

(gravitational, centrifugal, Coriolis, electromagnetic etc.).

For Newtonian fluids, the molecular rate of momentum transport can be written as:

T=−p+2

3µ∇ · uI+ 2 µD,(2.18)

where Tis the stress tensor, pis the static pressure, µis the dynamic viscosity, Iis the

unit vector and Dis the rate of deformation tensor. The rate of deformation tensor D

is defined as:

D=1

2h∇u+ (∇u)Ti.(2.19)

16 CHAPTER 2. MATHEMATICAL FORMULATION

The above two [Eq. 2.18 and Eq. 2.19] equations in tensor notation gives:

Tij =− p+2

3µ∂uj

∂xj!δij + 2 µ Dij ,(2.20)

Dij =1

2 ∂ui

∂xj

+∂uj

∂xi!,(2.21)

where δij is Kronecker symbol (δij =1 if i=jand δij=0 otherwise).

Finally, the integral form of the momentum conservation equation with the body forces

represented by fbecomes:

∂

∂t (ρu)dV +Z

ρu u ·ndS =Z

T·ndS +Z

ρf.(2.22)

By applying Gauss divergence theorem to the convective (2nd term) and diffusive (3rd

term) flux terms in the above equation, results:

∂(ρu)

∂t +∇ · (ρu u) = ∇ · T+ρf.(2.23)

The above equation represents a coordinate-free form of the momentum conservation

equation. The momentum conservation equation [Eq. 2.23] together with the conti-

nuity equation [Eq. 2.14] are popularly known as Navier-Stokes equations.

2.2.4 Species Equation

In addition to the continuity and momentum equations, if a problem involves transport

of mass then an additional species concentration equation needs to be solved. The

conservation equation for species transport takes the general form:

∂

∂t (ρ c)dV +Z

ρ c u·ndS =−Z

J·ndS , (2.24)

in the absence of a chemical reaction. Here, cis the specie concentration and Jis the

diffusion flux of the species. Diffusion flux arises due to the concentration gradients

and is defined by Fick’s first law of diffusion. According to this:

J=−D∂c

∂x ,(2.25)

2.2. GOVERNING EQUATIONS OF FLUID FLOW 17

where Dis the diffusion coefficient.

Present thesis do not deal with the problems related to heat transfer and hence the

energy conservation equation is not discussed here. It is clear from the above equations

that they all have common terms and it is useful to write the conservation equations

in one general form. Thus, the integral form of the generic transport equation with φ

as a variable gives:

∂

∂t (ρ φ)dV +Z

ρ φ u·ndS =−Z

Γ ∆ φ·ndS +Z

SφdV . (2.26)

The coordinate-free form of the above equation yields:

∂(ρ φ)

∂t +∇ · (ρ φ u) = −∇ · (Γ ∆ φ) + Sφ.(2.27)

The Eq. 2.27 highlights the various transport processes: the rate of change and the

convective terms on the left hand side whereas the diffusive and the source terms on

the right hand side. The next chapter deals with the numerical methods which are

used to solve the transport equations derived in this chapter.

18 CHAPTER 2. MATHEMATICAL FORMULATION

Chapter 3

Numerical Methodology

3.1 Introduction

Experimental methods are popular and desirable in science and engineering and they

also provide a deep understanding of important concepts about a particular problem.

But sometimes it becomes very difficult to establish an experimental model. Moreover,

the flows and related phenomena described by integro-differential equations cannot be

solved analytically except in some special cases. In such instances, numerical methods

are mandatory. Solution of transport equations using numerical method is known as

computation. The advantages of computation are low cost, faster than experimental

method, gives local values compared to experiments where only global values are

available and ability to simulate both ideal and real conditions. However, it might be

difficult to obtain numerical solution for problems involving complex geometries, non-

linearities, sensitive variations. In addition, it can be difficult to develop mathematical

model for some processes like combustion, multi-phase flows, non-Newtonian flows,

complex turbulent flows etc. Even if the problem is modelled, it is difficult to determine

how reliable is the modelled equation. Thus verification of the model has to be done

from the results of corresponding experiments. The fundamental equations of the fluid

flow have been derived in Chapter 2 and simplification of these equations are used to

realize the fluid flow. The simplified equations are usually based on a combination of

20 CHAPTER 3. NUMERICAL METHODOLOGY

approximations and dimensional analysis.

Inspite of the limitations, some numerical methods are used to obtain an approx-

imate solution for the transport equations. For this puprose, one need to follow a

discretization method which approximates the differential equations by a system of

algebraic equations which can be then solved on a computer. This method is called

Computational Fluid Dynamics (CFD). The first step in CFD is the development of

mathematical model which contains set of differential equations and boundary con-

ditions based on the conservation principles. Next step is the division of solution

domain into a small finite number of sub domains (numerical grid or descrete points)

at which the variables are calculated. This numerical grid is generated depending

on the problem and can be of structured or unstructured. A suitable discretization

method (finite difference, finite volume or finite element etc) is chosen to solve the

mathematical model at the descrete points in space and time. Discretization methods

applied at descrete points yields a system of non-linear algebraic equations and an

iterative method is used to solve these non-linear algebraic equations. Based on the

accuracy and efficiency of the solution, one needs to set the convergence criteria for

the iterative method. Finally, post processing of the solution into a physically realis-

tic result completes the numerical investigation of any fluid flow. Since discretization

plays very important role in numerical methodology, the following section is dedicated

to finite volume descitization method which is used in the present thesis.

3.2 Finite Volume Method

Basically, there are three well-known discretization schemes which are finite difference

(FD), finite volume (FV) and finite element (FE) available to solve a mathematical

model. Finite difference method is the oldest method available and the partial deriva-

tives are approximated by the nodal values of variables using Taylor series expansion

at each grid point resulting in algebraic equation. Finite difference method is very

simple and effective for simple geometries. In finite volume method whole domain is

3.2. FINITE VOLUME METHOD 21

discritized into small volumes and differential equation is integrated and terms are

approximated to obtain algebraic equation at each control volume. In finite element

method complete domain is divided into finite elements that are generally unstruc-

tured. In this method differential equation is multiplied by a weight function before

integration. Unknown variable is approximated by a shape function and substituted

in the weighted integral equation.

Finite volume method is popular among the engineers and scientists because each

term in this method has a physical meaning. The solution domain is divided into a

finite number of small control volumes called computational cells where the flow vari-

ables are defined. This is usually done at the geometric centre of the a control-volume

and this approach is called colocated cell-centered arrangement. The differential form

of the governing equations are first integrated over the individual computational cells

to transform the flux terms into surface integrals and then approximated in terms of

the cell-centred nodal values of the dependent variables.

The integral form of the generic transport equation with φas a variable is given

as:

∂

∂t (ρ φ)dV +Z

ρ φ u·ndS =Z

Γ ∆ φ·ndS +Z

SφdV . (3.1)

Eq. 3.1 applies to each control volume as well as to the complete solution domain. To

obtain an algebraic equation for each control volume, the surface and volume integrals

need to be approximated using discretization techniques. The discretization of each

term in the Eq. 3.1 is performed based on the name convention depicted in Fig. 3.1.

In the current formulation, each control volume centre has six immediate neighbours

which get involved in each term of the generic transport equation. The control volume

surface is subdivided into six plane faces denoted by e,w,n,s,t and b to represent east,

west, north, south, top and bottom directions with respect to the central node Pas

shown in Fig. 3.1.

22 CHAPTER 3. NUMERICAL METHODOLOGY

•P••

•

•B

•

•n

•s

•P••

•

•B

•

•n

•s

Figure 3.1: Illustration of control volume approach with three dimensional discretiza-

tion notation.

3.2.1 Spatial Discretization

The second and third terms of the Eq. 3.1 are contributions due to the convection

(ρ φ u·n)and diffusion (Γ ∆ φ·n). The net flux through the control volume boundary

is the sum of integrals over the six control volume faces, this means:

f dS =X

CSk

f dS , (3.2)

where fis the component of the convective or diffusive vector in the direction normal

to the control volume face. One would need to know fto calculate the surface integral

in Eq. 3.2 and therefore the cell face values are approximated in terms of the nodal

(cell centred) values. Similarly the source term in Eq. 3.1 requires integration over

the volume of a control volume. A second-order approximation is used to replace the

volume integral by the product of the mean and the volume of a control volume. Thus,

the volume integral is approximated as the value at the centre of a control volume as:

SP=Z

S dV =S∆V≈SP∆V , (3.3)

3.2. FINITE VOLUME METHOD 23

where SPis for the value of Sat the control volume centre. No interpolation is

necessary to calculate SPsince all variables are available at node P.

The approximations to the integrals also require the values of the variables at

locations other than control volume centres. Assume that u,ρand Γare known at

all locations. The value of φand its gradient normal to the cell face at one or more

locations on the control volume surface are needed to calculate convective and diffusive

fluxes. The volume integrals of the source term also require these values. Hence, φ

has to be expressed in terms of the nodal values by interpolation methods. There are

several interpolation methods availabal for approximating φand some of the frequently

used among them are described below.

Upwind Differencing Scheme (UDS)

In this interpolation method, φeis approximated using a forward or backward differ-

ence approximation depending on the flow direction. If the flow is in positive direction

(u·n)e>0, approximation for φegives:

φe=φP,

if the flow is in negative direction (u·n)e<0, approximation for φegives:

φe=φE.

For the flow in positive direction, Taylor series expansion about Pis given as:

φe=φP+ (xe−xP) ∂φ

∂x!P

+(xe−xP)2

2 ∂2φ

∂x2!P

+H , (3.4)

where Hdenotes higher-order terms. This approximation retains only the first term

on the right hand side of the Eq. 3.4, so it is a first order scheme. Although UDS ap-

proximation never yields oscillatory solutions, it is numerically diffusive due to Taylor

series truncation error.

Central Differencing Scheme (CDS)

Central differencing approximation for the value at control volume face centre is linear

interpolation between two nearest neighbouring nodes, irrespective of flow direction.

24 CHAPTER 3. NUMERICAL METHODOLOGY

Approximation for φeat location eis given as:

φe=φEλe+φP(1 −λe),(3.5)

where the linear interpolation factor λeis defined as:

λe=xe−xp

xE−xP

The Taylor series expansion for CDS gives:

φe=φEλe+φP(1 −λe)−(xe−xP)(xE−xe)

2 ∂2φ

∂x2!P

+H . (3.6)

This is the simplest second order accurate and is most widely used scheme. This

scheme may produce oscillatory solution because of the large variations of the fluxes

between the neighbouring points, but produces less numerical diffusion.

Quadratic Upwind Interpolation for Convective Kinematics (QUICK)

The quadratic upstream interpolation for convective kinetics (QUICK) scheme [41] is

a third order scheme which fits parabola through two points upstream and one point

downstream to get an interpolated value. Approximation for φeis given according to

the nature of convection. If the flow is in positive direction:

φe=φP+g1(φE−φP) + g2(φP−φW),

and if the flow is in negative direction:

φe=φE+g3(φP−φE) + g4(φE−φEE),

where the coefficients g1,g2,g3and g4are nodal coordinates. Details of the expressions

of these coefficients can be found in [21]. The Taylor expansion of this scheme gives:

φe=6

8φP+3

8φE−1

8φW−3(∆x)3

48 ∂3φ

∂x3!P

+H . (3.7)

The first three terms on the right hand side represent the QUICK approximation,

while the last term is the principal truncation error. As we can see, this quadratic

3.2. FINITE VOLUME METHOD 25

interpolation has a third order truncation error.

Diffusive term is generally discretized using CDS approximation whereas convective

term can be discretized by any scheme depending on the strength of the convection.

A dimensionless number, the Peclet number:

Pe =ρu

(Γ/δx),(3.8)

which is the ratio of strength of convection to diffusion, is introduced. At high Pe

number, results produced from CDS are oscillatory because convection is high and

the complete contribution comes from the upstream point, nothing comes from the

downstream point. Therefore if the flow is convection dominated then the upwind

differencing (backward or forward) scheme should be used. Note that QUICK scheme

should not be used for all-tetrahedral meshes and that it can be dispersive.

Finally, the spatial discrete equation takes the general linearised form:

aPφP=X

anbφnb +b . (3.9)

where aPdenotes the contribution of fluxes and other terms at the control volume

centre and anb are the coefficients of the neighbouring nodes (east, west, north, south,

top and bottom). The term bdenotes the source term of the discrete equation.

3.2.2 Temporal Discretization

Essentially, there are two temporal discretization schemes available to perform the

transient simulations and they are the first order fully-implicit scheme and second

order Crank-Nicolson scheme. Under fully-implicit formulation, the fluxes prevailing

over the time interval are calculated from the new time-level values of the variables.

The first term in Eq. 3.1 represents the time derivative and on discretization of it using

fully-implicit scheme gives:

T1 = (ρ φV )n+1

P−(ρ φV )n

∆t,(3.10)

26 CHAPTER 3. NUMERICAL METHODOLOGY

where the superscripts nand n+ 1 refer to old and new time levels respectively,

seperated by an interval ∆t.

In principle, the fully implicit scheme allows any magnitude of time step to be used, but

for transient problems, ∆tmust be small enough to limit the temporal approximation

errors to acceptable levels. In Crank-Nicolson scheme, the point at n+ 1/2is used as

the expansion point for the second-order formulation. This scheme is liable to generate

unbound extrema and this trend can be improved by the use of smaller ∆t.

3.2.3 Boundary Conditions

All CFD problems are defined in terms of initial and boundary conditions. In order to

solve the set of algebraic-equations, one need to estimate the boundary fluxes either

explicitly or in terms of the variables at inner nodes. The most commonly used

boundary conditions have been described here.

•Inlet: The realization of this boundary condition is relatively easier, where the

destribution of mass flux and fluid properties are known.

•Outlet: The total mass flux from all inlet planes should be equal to the total

mass flux from all outlet planes. This means, the gradients of all variables along

the flow direction at the outflow surface are taken to be zero and the exit mass

flow is fixed from overall mass conservation criterion.

•Symmetry: The symmetry or zero-gradient boundary condition is a common

type of boundary condition in all CFD solvers. The conditions at the symmetry

boundary are, no flow and no scalar flux across the boundary. In the implemen-

tation, normal velocities and normal gradients of all other variables are set to

zero at a symmetry boundary.

•Slip surface: A slip surface is impermeable to flow across it but offers no resis-

tance to tangential motion and is usually appropriate for inviscid flow calcula-

tions and it can also be used for free surface flows. Slip condition is implemented

3.3. ROTATING AND MOVING MESHES 27

by setting the mass and tangential momentum fluxes to zero at cell faces as well

as zero normal velocity and zero normal gradients for all other variables except

pressure.

•Wall: The wall boundary condition is the most encountered boundary condition

in general fluid flow problems. The no-slip boundary condition is applied directly

in the case of laminar flows and certian turbulence models via wall functions is

used in the case of turbulent flows. A wall moving at a known velocity may also

be prescribed to account for rotating flow problems.

3.3 Rotating and Moving Meshes

In the present thesis as well as in many practical applications of fluid dynamics, fluid

motion is induced or governed by relative movement between stator and rotor. This

is usually accompanied by the strong characteristic unsteadiness in the flow pattern.

Important examples where such flows occur are:

• kneaders, single and twin screw extruders

• mixing vessels

• rotary and reciprocating engines

• turbomachinary, axial and centrifugal pumps

• ship and aircraft propellers.

The mostly used mixing equipment is stirrer. In the case of axi-symmetric geometries,

the computational domain can be made fixed in space using rotating coordinate sys-

tem. The similar idea is applied in case of sliding interface technique for non-symmetric

vessels. In this case the computational domain is divided into partial domains and the

domains with the rotating parts are calculated using rotating reference frame. Several

CFD codes offer modeling of flows involving rotating bodies. In the present work, a

28 CHAPTER 3. NUMERICAL METHODOLOGY

finite volume based commercial code Star-CD from CD-adapco is used to investigate

flows in stirred vessels. Major features available in Star-CD for rotating geometries

are rotating reference frame method and general mesh motion [1]. A brief discussion

of these features is given in this section with an example.

3.3.1 Rotating Reference Frame

In the problems of stirred vessels, the impeller is generally rotated with certain speed.

The influence of rotation can be realized either by prescribing tangential speed at

the boundary cells or by transforming the equations to a rotating frame of reference.

Velocities and accelerations recognized by an observer in a rotating frame of reference

are different from that recognized by an observer at rest in an inertial or fixed frame of

reference. To futher elaborate, consider a fluid in rotation about some fixed axis with

an angular velocity (Ω). Let (x1,y1,z1) be the fixed coordinate system and (x2,y2,z2)

be the rotating frame of reference. The fixed axis of rotation is common to both

the coordinate systems, i.e.,z1and z2coincide. The transformations relating the two

coordinate systems are:

x1=x2cos(Ωt) + y2sin(Ωt)

y1=y2cos(Ωt)−x2sin(Ωt)

z1=z2,

and,

x2=−x1sin(Ωt) + y1cos(Ωt)

y2=−y1sin(Ωt) + x1cos(Ωt)

z2=z1.

In order to calculate velocities and accelerations in rotating reference frame, again

consider a fluid element with position vector rin inertial frame of reference and observe

the motion of this fluid element in rotating reference frame which rotates with a

3.3. ROTATING AND MOVING MESHES 29

constant angular velocity (Ω). Then, the equation of motion in inertial frame is

obtained by time derivative on rand read as:

dt =Ω×r.

The time derivative of rin rotating reference frame has two components, one from

inertial reference frame and another from its own rotation. The above equation is

re-written in rotating reference frame by adding instantaneous velocity u0as:

dt =u=u0+Ω×r.

The second time derivative operating on position vector rgives an acceleration as:

a= d2r

dt2!= du

dt !

= d

dt0+Ω×!(Ω×r)

After simplification and re-arranging the above equation yields:

a0=a−Ω×(Ω×r)

|{z }

−2Ω×u0

|{z }

,(3.11)

where the terms I=Ω×(Ω×r)and II = 2Ω×uare centrifugal and Coriolis

accelerations respectively. The centrifugal and Coriolis accelerations result in fictitious

forces in rotating reference frame. Adopting the above mentioned transformation

procedure, the momentum equations in rotating reference frame are modified by taking

into account Coriolis and centrifugal force terms. The Navier-Stokes equations written

for use in a rotating reference frame at constant angular velocity (Ω) then read as [80]:

∂u

∂t + (u· ∇)u+Ω×(Ω×r) + 2Ω×u=−1

ρ∇p+ν∇2u+1

ρb.(3.12)

Note that in the above Eq. 3.12 the prime symbol 0is deliberately omitted. The

continuity equaion (∇ · u= 0) is invariant under the given transformation.

The rotating reference frame feature of Star-CD enables to model cases where the

entire mesh is rotating at a contant angular velocity around a prescribed axis as well as

different mesh blocks moving with different angular velocities around different rotating

30 CHAPTER 3. NUMERICAL METHODOLOGY

axes. This means, there are two types of rotational motion that may be simulated,

corresponding to:

• Global rotation, using single rotating reference frame.

• Region wise mesh rotation using multiple reference frames.

Global Rotation

This is the simple rotation where the entire mesh rotates with a prescribed angular

velocity (Ω) about a given axis. The main features of this type are:

• The axis of orientation and sense of rotation are defined by the z-axis of a local

coordinate system and the convention that Ωobeys right hand rule.

• Selection of this feature adds the Coriolis and centrifugal force terms to the

momentum equations and yields Eq. 3.12.

Region Wise Mesh Rotation

This feature allows different mesh blocks within the model to rotate at different angular

velocities, i.e. the rotational source terms added to the relevant equations depending

on the local spin and axis of rotation. The basic idea here is to calculate relative mesh

rotation by activating the appropriate rotational terms at cells within each specified

region, without rotating the entire mesh. Then it is possible to have multiple adjoining

rotating regions, each with its own rotational speed.

3.3.2 Mesh Motion

In a number of situations it is convenient for the computatonal grid to follow moving

boundaries rather than rotating the entire mesh. In this case of moving grids, the

time coordinate has to be transformed in a way similar to the space transformation.

A brief description of a moving grid system is given below. On the basis of the rules of

3.3. ROTATING AND MOVING MESHES 31

partial differentiation, the total time derivative from the Lagrangian view point reads:

dφ

dt =∂φ

∂t +∂φ

∂xi

dxi

dt ,(3.13)

where the left hand side is the time rate of change of φof a given fluid element when it

moves from one point to the other point. The first term on the right hand side is the

time rate of change of φat a fixed point. The quantity ∂φ

∂xiis the time rate of change

of φat a given instant of time and dxi

dt is the grid velocity ug. Therefore, Eq. 3.13 can

now be written as:

∂φ

∂t =dφ

dt −ug

∂φ

∂xi

.(3.14)

The above Eq. 3.14 can be referred as Arbitrary Lagrangian-Eulerian (ALE) formula-

tion for moving grids. The conservation form of the ALE is obtained with the help of

the geometric or space conservation law (SCL), which ensures that the rate of change

of volume is consistent with net volume flux due to grid movement [18]:

dt Z

φ dV −Z

φub·ndS = 0 .(3.15)

The generic transport Eq. 3.1 in finite-volume formulation using Eq. 3.15 for moving

grids yield:

dt Z

ρ φ dV +Z

ρ φ (u−ub)·ndS =Z

Γ ∆ φ·ndS +Z

SφdV . (3.16)

The above equations are used if the control volume does not move. If the boundary

moves with the same velocity as the fluid, the mass flux through the control volume

face will be zero. If this is true for all control volume faces, then the same fluid remains

in the control volume and becomes a control mass. The resulting conservation equation

describes the Lagrangian fluid motion.

Star-CD offers several methods to carry out moving mesh simulations. Of course,

the mesh motion is not entirely arbitrary; there are certain limits on the degree of

distortion that can be tolerated, imposed on accuracy and stability requirements. The

moving mesh feature in Star-CD is activated by command MVGRID and changes in

32 CHAPTER 3. NUMERICAL METHODOLOGY

grid movement can be specified either by change grid operation in the EVENTS com-

mand module, or by user defined function included in subroutine NEWXYZ. In this

subroutine, the geometry of the model can be varied by defining vertex coordinates as

the function of time. There are several methods to apply above mentioned operations.

They are:

• Cell layer addition or removal.

• Sliding mesh.

• Cell attachement and change of fluid.

• Mesh region inclusion or exclusion.

Sliding mesh method is utilized in the present thesis and therefore, it is described here.

The detailed information of other methods can be found in [1].

Sliding Mesh

Various methods of implementing sliding mesh are regular sliding interface, arbitrary

sliding interface (ASI) and automatic events generation. The regular sliding inter-

face method combines both the cell attachment and the change grid operation of the

EVENTS command module. The ASI approach involves re-computing boundary face

matches at each time step, and may thus be slightly slower than using ATTACH

events. Automatic events generation is a simplified procedure for generating a tran-

sient moving mesh model, starting from a stationary mesh built at a given position of

the impeller blade. The events, moving grid commands, attach boundaries and any

other required user subroutines are generated in a single command. The next section

expands the rotating reference frame and sliding (moving) method with the aid of a

simple example.

3.3. ROTATING AND MOVING MESHES 33

(a) Vertical (b) Horizontal

Figure 3.2: Schematic of the simple reference geometry considered for rotating and

sliding mesh analysis.

3.3.3 Simple Reference Geometry

The purpose of using simple reference geometry is to apply and compare rotating

reference frame simulations with moving mesh simulations. The geometrical configu-

ration used in the present study contains two concentric cylinders height 7 m, outer

cylinder radius 3 m and inner cylinder radius 0.5 m. Each side of the cube is 2.5 m and

is attached to the inner cylinder (shaft) as shown in Fig. 3.2.

Simulation

The simulation methods employed for the simple reference geometry are the rotating

reference frame and the moving mesh with change grid operation. In the first case, the

entire mesh is rotated with an angular velocity of 10 rpm. Since the boundaries of the

rotating domain assumed to be rotating with the same angular velocity, an opposite

spin of the same angular velocity is defined on inner walls to make them stationary. In

the moving mesh case where the grid (cell vertices) is changed at each time step using

change grid operation with events. No-slip boundary conditions are used at the outer,

34 CHAPTER 3. NUMERICAL METHODOLOGY

top, bottom and inner walls. All fluid motion is caused by the rotation of the cube and

shaft. The unstructured mesh contains approximately 50000 grid nodes is generated

using Star-Design and the numerical simulations are performed with Star-CD v 3.22.

A time step of 0.01s is used in the simulations and 1200 time steps were performed in

each case corresponding to two rotations of the shaft. Estimated calculation time is

one hour per rotation on an AMD Athlon XP, 3200+ and 1GB RAM single processor.

(a) After 6 sec (one rotation) (b) After 6 sec (one rotation)

Figure 3.3: Section plots of velocity vectors after one rotation of the impeller in moving

mesh (left) and rotating reference frame (right) simulations.

Results and Discussions

Fig. 3.3 shows the section plots of velocity magnitudes after one rotation of the impeller

in the simple reference geometry. The flow pattern is shown by means of velocity

vectors. The length of the vectors is proportional to the magnitude of the liquid

velocity. As expected, two circulation loops can be seen around the impeller in both

cases.

Fig. 3.4 shows the section plots of the velocity magnitude fields in both moving

3.3. ROTATING AND MOVING MESHES 35

(a) After 3 sec (b) After 3 sec

(e) After 12 sec (f) After 12 sec

Figure 3.4: Contour section plots of velocity magnitudes at different time intervals in

both moving (left) and rotating reference frame (right) simulations.

36 CHAPTER 3. NUMERICAL METHODOLOGY

mesh and rotating reference frame simulations at different time steps. The flow pattern

is shown by means of contour plots. Since the motion of the fluid is governed by the

movement of the impeller (cube and shaft together), the highest velocities are found

near the impeller tip and the velocities are lower near the outer cylinder wall. Moving

mesh simulation results are shown on the left side while the rotating reference frame

simulation results are shown on right side. In Fig. 3.4 (a) and (b) are observed after 3

sec, i.e. after half of the roation of the impeller. As one can see, the flow is very week

away from the impeller. The velocity profiles are shown after one rotation in Fig. 3.4 (c)

Figure 3.5: Computational sensor points to extract quantitative information for better

comparison of the velocity magnitudes.

and (d). At this instance, the flow field becomes stronger and developing towards the

wall of the outer cylinder. After two rotations, much stronger flow fields are observed

in Fig. 3.4 (e) and (f). The principle reflection in both the moving mesh and rotating

reference frame simulations is that the results show very good coincidence with each

other. For more quantitative validation of the simulation results, 20 computational

sensing points are considered starting from the impeller wall to the outer cylinder wall

in radial direction as shown in Fig. 3.5. A very good agreement is noticed for for

3.4. TORQUE COMPUTATION IN A COUETTE-FLOW 37

0.4

0.8

1.2

1.6

0.0 0.4 0.8 1.2 1.6 2.0

Distance (m)

Velocity magnitude (m/s)

Mov-3sec

Rot-3sec

Mov-6sec

Rot-6Sec

Mov-12sec

Rot-12Sec

Figure 3.6: Quantitative comparison of velocity magnitudes computed both in moving

mesh and rotating reference frame simulations.

the simulation results of both the moving mesh and the rotating reference frame as

depicted in Fig. 3.6. It can also be seen that as the time proceeds, more uniform profiles

are achieved. From the above simulations, it is evident that both the moving mesh

and the rotating reference frame simulation give similar results. However, the moving

mesh calculations needs longer computational time than the rotational reference frame.

Therefore, the rotating reference frame method is utilized to accurately predict the

time dependant flows in stirred vessels. The following section explores the utilisation

of Star-CD to compute torque in rotating system by considering a simple Couette-flow

between two concentric rotating cylinders.

3.4 Torque Computation in a Couette-Flow

3.4.1 Introduction

Couette-flow refers to the laminar flow of a fluid between two concentric cylinders, one

of which is moving relative to the other. Couette-flow reactors are the widely used

38 CHAPTER 3. NUMERICAL METHODOLOGY

devices in chemical process industries especially due to their mixing properties. In this

section, the torque computation in a Couette-flow reactor using Star-CD are presented.

These computations are compared with the theoretically calculated torque as well as

torque computed using other well known commercial code Fluent. The simulations

and calculations are carried out for the different rotational speeds of the outer cylinder.

In the present investigations, a steady tangential laminar flow is considered between

(a) (b)

Figure 3.7: Schematic of (a) the Couette-flow geometry and (b) the computational

grid which contains 207,000 hexahedral control volumes.

two vertical concentric cylinders shown in Fig. 3.7. The fluid with constant density

and viscosity is used as a working fluid. The analytical solution for torque when outer

cylinder rotates with an angular speed Ω, is given as [22]

τ= 4πµLΩr2

2"(r2

1/r2

1−(r2

1/r2

2)#.(3.17)

Where τis the torque in N-m, µis the viscosity of the liquid in kg/m-s, Lis the height

of the cylinders in m, r1is the radius of the inner cylinder in m and r2is the radius

of the outer cylinder in m. The diameters of the inner and outer cylinders are 71 mm,

80.7 mm respectively and both have the same height of 9.7 mm.

3.4. TORQUE COMPUTATION IN A COUETTE-FLOW 39

Table 3.1: Comparison of torques calculated theoretically and using Fluent.

Angular Torque Torque Percentage

V elocity (Theoretical) (Fluent)Deviation

rad/s N −m N −m%

0.1 6.82E-08 6.82E-08 0.00

0.2 1.36E-07 1.37E-07 0.73

0.3 2.05E-07 2.05E-07 0,00

0.4 2.73E-07 2.74E-07 0.37

0.5 3.41E-07 3.43E-07 0.59

0.55 3.75E-07 3.78E-07 0.79

0.6 4.09E-07 4.13E-07 0.97

0.65 4.43E-07 4.48E-07 1.12

0.E+00

1.E-07

2.E-07

3.E-07

4.E-07

5.E-07

0 0.15 0.3 0.45 0.6 0.75

Angular Velocity (rad/s)

Torque (N-m)

Theoretical

Fluent

Figure 3.8: Comparison of theoretically calculated torque with torque computed using

Fluent.

40 CHAPTER 3. NUMERICAL METHODOLOGY

Table 3.2: Comparison of torques calculated theoretically and using Star-CD.

Angular Torque Torque Percentage

V elocity (Theoretical) (Star −CD)Deviation

rad/s N −m N −m%

0.1 6.80E-08 5.82E-08 14.4

0.2 1.36E-07 1.26E-07 7.35

0.3 2.04E-07 1.94E-07 4.9

0.4 2.72E-07 2.57E-07 5.51

0.5 3.40E-07 3.30E-07 2.94

0.55 3.74E-07 3.67E-07 1.87

0.6 4.08E-07 3.99E-07 2.21

0.65 4.42E-07 4.31E-07 2.49

0.E+00

1.E-07

2.E-07

3.E-07

4.E-07

5.E-07

0 0.15 0.3 0.45 0.6 0.75

Angular Velocity (rad/s)

Torque (N-m)

Theoretical

Star-CD

Figure 3.9: Comparison of theoretically calculated torque with torque computed using

Star-CD.

3.4. TORQUE COMPUTATION IN A COUETTE-FLOW 41

Table 3.3: Comparison of torques calculated theoretically and using Star-CD for highly

viscous liquid.

Angular Torque Torque Percentage

V elocity (Theoretical) (Star −CD)Deviation

rad/s N −m N −m%

0.3 2.04E-03 2.05E-03 0.49

0.35 2.38E-03 2.39E-03 0.42

0.4 2.72E-03 2.73E-03 0.37

0.45 3.06E-03 3.07E-03 0.33

0.5 3.40E-03 3.42E-03 0.58

0.55 3.74E-03 3.76E-03 0.53

0.6 4.08E-03 4.10E-03 0.49

0.65 4.42E-03 4.44E-03 0.45

0.E+00

1.E-03

2.E-03

3.E-03

4.E-03

5.E-03

0 0.15 0.3 0.45 0.6 0.75

Angular Velocity (rad/s)

Torque (N-m)

Theoretical

Star-CD

Figure 3.10: Comparison of theoretically calculated torque with torque computed

using Star-CD for highly viscous liquid.

42 CHAPTER 3. NUMERICAL METHODOLOGY

3.4.2 Simulation

To compute the torque in couette-flow reactor, numerical calculations were performed

using the two widely used commercial CFD softwares Star-CD and Fluent. These

computations are compared with the theoretically calculated values. A viscosity of

0.001 kg/m s is used in the calculations of Star-CD and 0.00103 kg/m s is used in the

case for Fluent calculations.

Star-CD

A 3-dimensional block structured grid with 207,000 computational cells were gener-

ated using Pro-STAR. A steady, laminar, isothermal simulations were performed at

different Reynolds numbers corresponding to different angular velocities. The torque

was computed once the convergence is obtained. The detailed method of calculating

torque using Star-CD is given in section 4.2.2.

Fluent

For the simulations using Fluent, the results were taken from [29]. Wall shear stress

over the area is calculated after convergence is achieved. Integration of wall shear

stress over the area of the inner cylinder and multiplied with the the radius of the

inner cylinder gives the torque required to rotate the inner cylinder.

3.4.3 Results and Discussions

The theoretical values for the torque and the torque calculated using Star-CD are

tabulated in Table 3.2. Fig. 3.9 shows a comparison between the theoretical value

and the torque calculated using Star-CD. It can be observed in this figure that there

exists a small discrepancy in these values. Where as the variation between the torque

[Table 3.1, Fig. 3.8] obtained using Fluent and the theoretically calculated torque are

in good agreement. After a carefull observation of the results, it can be seen that the

percentage of deviation of the torque values obtained using Fluent from the theoretical

3.4. TORQUE COMPUTATION IN A COUETTE-FLOW 43

torque values is in the order of one. But the percentage of deviation of the torque

values calulated using Star-CD from the theoretical torque values is in the range 2-

15. From this we can be seen that the Fluent is more appropriate than Star-CD for

this particular case of the torque computation in a couette-flow reactor. However,

the intensity of deviation decreases when highly viscous liquids are considered in the

investigations. These findings can be found in Table 3.3 and Fig. 3.10. It is believed

from these observations that Star-CD predicts values close to theoretical values at

high viscosity of the liquids.

44 CHAPTER 3. NUMERICAL METHODOLOGY

Chapter 4

CFD-Analysis of Anchor Mixer

4.1 Introduction

Fluid mixing in a stirred vessel agitated with an anchor impeller is discussed in this

chapter. The anchor impeller has a simple and basic configuration suited for mixing

of highly viscous liquids, and is widely used in the chemical and food industries. The

principle advantages of anchor mixer are that it prevents sticking of pasty materials,

and promotes good heat transfer with the wall because of close clearance between the

impeller and vessel wall. Because of its unique advantages, the detailed analysis of the

flow around the vertical blades of an anchor impeller gains much appreciation. This

analysis includes information regarding power consumption, mixing and heat transfer.

Some of the early experimental studies on flows in an anchor mixer were reported

in [9,60]. They provided the detailed description of the power input and flow fields with

Newtonian and non-Newtonian fluids. Early two-dimensional numerical investigations

were developed by [39] using an iterative method. This paper reports the flow of

highly viscous fluids in anchor mixer and gives information on flow characteristics.

Extensive three-dimensional numerical analyses were then followed during early 90s

[2, 10, 33, 34, 59, 74]. Inspite of the importance of anchor mixer in process industries,

there are only few studies in the literature which deal with this mixer type. The

majority of works for stirred vessels refer to six-blade Rushton turbines and paddle

46 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

Figure 4.1: Schematic representation of anchor mixer, where all the dimensions are in

mm (Geometrical configuration is received from BASF AG, Ludwigshafen).

type impeller mixers under turbulent flow condition. This work contributes for the

simulation of anchor mixer.

4.2 Analysis of Partly Filled Anchor Mixer

The present section deals with a three-dimensional analysis of a partly filled anchor

mixer which is shown in Fig. 4.1. Since the vessel is partly filled, Volume of Fluid

(VOF) method is used to account for the free surface flow. One of the leading com-

mercial CFD codes Star-CD from CD-Adapco has been used for flow field analysis,

torque computation and power number calculations. The flow considered here is lam-

inar, Newtonian and time dependant in a two arm anchor type stirrer fixed to a shaft

and placed at the centre of the vessel. Simulation procedure addresses the numerical

solution of the three-dimensional Navier-Stokes equations in rotating reference frame

coordinates with additional centrifugal and Coriolis force terms [Eq. 3.12].

4.2. ANALYSIS OF PARTLY FILLED ANCHOR MIXER 47

4.2.1 Volume of Fluid (VOF) Method

Since the mixer is partly filled with the liquid, the free surface flow has to be accounted

for the analysis of two-phase mixing inside the anchor mixer. Therefore, this problem

is solved using Volume of Fluid (VOF) method to implicitly track the residence region

of two immiscible flows. According to this, the interface between the liquid and the gas

is captured at successive time steps and is represented by the distribution of a passive

scalar αl, which is defined as the ratio of volume of liquid to the total volume of the

fluid (liquid and gas) in a computational cell. It was further assumed that the fluids

are incompressible and have a considerable density difference. Based on phase-related

conservation, advection of the disperse phase is governed by an additional transport

equation for the volume fraction αlof this phase:

∂αl

∂t +∇ · (αlu) = 0 ,(4.1)

determining the temporal propagation of the phase with respect to the flow field.

Hence, the balance equations of the mass and momentum in addition to Eq. 4.1 can

then be solved for both phases simultaneously. Using VOF method, the disperse liquid

phase corresponds to the region where αlhas the value 1 and the continuous phase

has the value of 0, while the interface is located within the grid cells for which 0 <

αl< 1. Therefore, cells with 0 < αl< 1 comprise both the disperse and continuous

phase. To further clarify this approach, Fig. 4.2 represents the distribution of αl

inside the calculation domain. In order to maintain the interface sharpness, Eq. 4.1 is

solved using a higher order compressive scheme, the so-called Compressive Interface

Capturing Scheme for Arbitrary Meshes (CICSAM). A requirement for this scheme is

that the Courant number Cu of the flow should not exceed a certain limit Cu∗(0.3

in this case). Otherwise the scheme is not able to maintain a sharp interface. The

Courant number Cu is defined as:

Cu =U·∆t

∆x.(4.2)

48 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

Figure 4.2: Cell occupation of liquid phase (0 < αl< 1) in a surrounding gas phase

(αl=0) based on the VOF-method.

The physical properties (density ρand viscosity µ) of the mixture are expressed as

functions of the phase properties and are as follows:

ρ=ρlαl+ρg(1 −αl),(4.3)

µ=µlαl+µg(1 −αl).(4.4)

4.2.2 Problem Specification and Numerical Approach

For the analysis, a concave-bottom cylindrical vessel with an anchor type stirre fixed to

a shaft and placed at the centre of the vessel is considered. An unstructured tetrahedral

mesh Fig. 4.3 is generated using ICEM-CFD and exported to Star-CD for numerical

analysis. The total number of computational cells used in these investigations are

approximately 300,000. The simulations were perfomed with different mesh refinement

and sufficient care has been taken in setting the solver parameters such as discretization

schemes, time step size etc. A viscosity of 27.25 Pas and density of 1120 kg/m3are

assigned for the liquid phase and air is considered as the gas phase with viscosity

1.8551E-05 Pas and density of 1.184 kg/m3. The vessel is filled initially with the

liquid upto a level of 148 mm and the stirrer rotates at an angular speed of 50 rpm.

4.2. ANALYSIS OF PARTLY FILLED ANCHOR MIXER 49

Figure 4.3: Computational grid which consists of approximately 300,000 tetrahedral

control volumes generated using ICEM-CFD tool. Mesh with magenta colour (below)

correspond to liquid phase and the one with aqua colour (above) correspond to gas

phase.

After solving the transport equations for velocity and pressure fields, torque and power

numbers are calculated based on the the formulae given in the next sections.

Torque

Torque or the moment of force is an axial vector quantity that represents the magnitude

of rotational or angular force applied to a rotational system at a distance from the axis

of rotation. This force is determined by linear force multiplied with the radius. More

generally, torque can also be defined as the rate of change of angular momentum.

Mathematically, the torque on a particle which has the position rin some reference

frame, can be represented as:

τ=r×Ft.(4.5)

Alternatively, if a force of magnitude Fis at an angle θfrom the distance rwithin

50 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

the plane perpendicular to the rotation axis, then from the definition of cross product,

the magnitude of the torque arising is:

τ=rF sin θ . (4.6)

The expression for computing flow induced forces on the wall cell faces are given

below. Shear force:

Fs=τwAbt,(4.7)

pressure force:

Fp=pbAbnb,(4.8)

total force:

Ft=Fs+Fp,(4.9)

and hence the total torque is given by:

τ=r×Ft,(4.10)

Power Number

The information provided in this section is based on [14]. The mechanical energy

balance is to be analysed for the first estimation of the mixing efficiency. It should be

noted that this mechanical energy is transferred in two different ways. On one hand,

the mechanical energy can be transferred through the inner boundaries, i.e., through

the agitation with the moving wall (ex. stirrer, extruder and kneader) and on the other

hand the mechanical energy can be transferred over outer boundaries i.e., through

pressure drops along the flow direction (ex. static mixer). Scalar multiplication of the

Navier-Stokes equations by the velocity, and integration over the considered control

volume V(t)in the mixing domain gives the following relation:

EKin =PF low +PAgit −PDiss .(4.11)

4.2. ANALYSIS OF PARTLY FILLED ANCHOR MIXER 51

The power transferred due to the flow:

PF low = ∆p˙

V+Z

Ain

2|u|2VaxdA −Z

Aout

2|u|2VaxdA , (4.12)

The power induced due to the mechanical agitation or mixing:

PAgit =Z

Amix(t) ηu·∂u

∂n−pu ·n!dA , (4.13)

and the dissipated power:

PDiss =ηZ

|∇u|2dV . (4.14)

Here, Amix(t)denotes the surface area of the moving inner boundary (mixing part),

∆pis the pressure drop between the inlet and outlet, ˙

Vis the volume flux through the

mixing apparatus and Vax is the component of velocity in axial direction. Eqs. 4.12-

4.14 are true under the following consideration that the inlet and outlet surfaces are

planar and their velocity vectors have to be perpendicular to them. Additionally, inlet

and outlets are governed by constant pressure. For simplification, the formulation is

considered only for Newtonian fluids. In the general case, shear stress due to the stress

tensor have to be considered [72].

Depending on the mixing process, either Eq. 4.12 or Eq. 4.13 have to be considered.

In the case where the moving walls exists, the considered volume domain depends on

time. In quasi steady state, where the kinetic energy is constant and hence Eq. 4.11

can be reduced to:

PDiss =PF low +PAgit .(4.15)

Further simplification of the Eq. 4.15 depends on the specific application. For batch

mixing processes in a mixing equipment with moving parts, PF low can be neglected.

Where as in the pressure driven flow mixing with agitation PAgit can be ignored. The

term on left hand side in Eq. 4.15 corresponds to the dissipation term and is a required

quantity for the integral estimation of mixing efficiency due to convection. For the

case, where the flow can be calculated, this dissipation term can be determined using

52 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

CFD-simulations. Depending on the complexity of the geometry, topology and scaling

of the mixing equipment, valid numerical calculation of the complete flow field can

be very demanding. Unfortunately, such kind of calculation even with the present

computational advancements cannot be possible. This is true for complex geometries

with self-cleaning and intermeshing (e.g. kneaders, extruders) design mechanisms.

Generally, computation of mixing mechanisms in extruders are difficult. However,

computation of single screw extruder can be done using appropriate moving or rotating

reference frame by considering as a channel flow [56]. Twin screw extruders cannot

be modelled using the sliding mesh technique because of the ovelapping of the two

screws. The fictitious domain method developed by [27] during the last years and

is shown appropriate [26]. The flow calculation in kneaders is even more difficult

because of the close clearance between the kneading pins. Recently, new techniques of

overlapping meshes have been developed [21, 30]. On the basis of Chimera methods,

individual meshes are devoted for each single moving objects. These moving meshes

are linked by a common back ground mesh. Because of above mentioned difficulties,

experimental data is important to develop correlations. As an integral parameter, the

total dissipated power PDiss is the fundamental parameter. In order to obtain a valid

correlation depending on the type of the mixing equipment, the parameters influencing

the power have to be determined. For example, in case of paddle impeller in which

the shearing term in Eq. 4.13 can be neglected and Eq. 4.13 can be simplified as:

ηZ

|u|2dV =Z

Amix(t)

pu·nmixdA ,(4.16)

where nmix is the outward normal vector of the mixing surface. For a thin pad-

dle impeller, the product of the pressure drop between the blades with the normal

component of the velocity has to be integrated. The pressure drop:

∆p=ρ

2U2∼ρ

2d2n2,(4.17)

with the blade diameter dand the rotational speed ncan be estimated. Since the

4.2. ANALYSIS OF PARTLY FILLED ANCHOR MIXER 53

normal velocity is proportional to dn and the surface is also proportional to d2, follows:

Amix(t)

pu·nmixdA ∼ρd5n3.(4.18)

This gives the following equation for the dissipated power:

PDiss =Neρd5n3,(4.19)

with the dimensionless number Ne called Newton number, relating the resistance force

to inertia force. According to Eq. 4.19, higher Newton number is advantageous because

the same power is dissipated by a smaller rotational speed.

If a force is allowed to act through a distance, it is doing mechanical work. Similarly,

if torque is allowed to act through a rotational distance, it is doing work. However,

time and rotational distance are related by the angular speed where each revolution

results in the circumference of the circle being travelled by the force that is generating

the torque. This means that torque that is causing the angular speed to increase is

doing work and the generated power may be calculated as:

P= 2πNτ,(4.20)

where Nis the rotational speed. From Eq. 4.19 and Eq. 4.20, we obtain the relation

between power number and torque as:

Ne =2πN

ρd5n3τ.(4.21)

4.2.3 Results and Discussions

For the preliminary investigations in the partly filled anchor mixer, calculations were

performed for two complete iterations. The results presented here correspond to the

time instant after one complete rotation of the anchor. Fig. 4.4 shows the free

surface flow field with computational grid and anchor stirrer at an instant of time.

A corresponding horizontal cross section with flow field and a rearrangement of the

interface is illustrated in Fig. 4.5. Here, it is clearly seen that the stirrer takes the

54 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

Figure 4.4: Computational grid together with the mixer geometry at the interface

between liquid and gas at an instant of time.

Figure 4.5: Computational grid at the interface between liquid and gas at an instant

of time.

4.2. ANALYSIS OF PARTLY FILLED ANCHOR MIXER 55

Figure 4.6: Representation of concentration of scalar variable in a mixer vessel stirred

with anchor.

interface along with it while moving there by the exchange between the two phases

makes possible. Fig. 4.6 demonstrates the contour plot of the volume fraction in

horizontal cross section. This profile shows a sharp interface between the liquid and

gas phase. Since the stirrer rotates counter clock-wise direction, one can notice the

crest on the front side and the trough behind the stirrer. Velocity contours along

the two-dimensional vertical and horizontal cross sections at the centre of the vessel

are shown in Fig. 4.7. These contours show a maximum velocity field around the

stirrer arms and a minimum velocity field outwards the centre of the vessel, which are

expected from the physical nature of the flow. These high velocity gradients results

in the flow field yields higher shear rates, which ultimately are responsible for reliable

mixing.

The pressure distribution inside the mixer is presented in Fig. 4.8. The high

pressure arises at the front of the stirrer and the low pressure behind the stirrer

and also in the rest of the regions inside the whole domain. Fig. 4.9 demonstrates the

evolution of torque with respect to time for a stirrer speed of 50 rpm. As time proceeds

at the initial state, torque decreases abruptly and reaches a constant value. At the

56 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

Figure 4.7: Contour plots of the velocity magnitudes in a mixed vessel stirred with

anchor. High velocity magnitudes can be seen arround the anchor arms.

4.2. ANALYSIS OF PARTLY FILLED ANCHOR MIXER 57

Figure 4.8: Pressure distribution in anchor mixer. High pressure results at the front

of the stirrer and low pressure behind the stirrer.

58 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

0.00

0.60

1.20

1.80

2.40

3.00

0.00 0.10 0.20 0.30 0.40 0.50

Time (s)

Torque (N-m)

Visc. 27.25 Pa.s

Figure 4.9: Time evolution of torque in anchor mixer.

onset of motion (t= 0), torque is at the maximum due to vessel resistance against

the fluid and stirrer motion. Torque reaches constant within a short time at high

viscosities, which is clearly observed in Fig. 4.10. Power number is expected to show

the same trend in its behaviour, since power or Newton number is the product of torque

and some constant quantity as shown in Eq. 4.21. There exists several ways for the

numerical investigation of mixing in stirred vessels. The well known among them are

numerical tracer experiments, Lagrangian particle tracking and entroy based measures.

Though numerical tracers experiments are widely used technique, it induces numerical

diffusion when dealing with highly viscous liquids. To avoid this, the present studies

were carried out based on Lagrangian particle tracking. In the preceeding section, the

highlights of this approach is emphasized and elaborated.

4.3 Particle Tracking

While numerical tracer experiments are of increasing importance as a means to analyze

mixing processes, there is a principle problem which becomes severe obstacle in case

of highly viscous liquid mixing or, more precisely, in case of highly schmidt numbers.

4.3. PARTICLE TRACKING 59

0.00

0.50

1.00

1.50

2.00

2.50

3.00

0.00 0.10 0.20 0.30 0.40 0.50

Time (s)

Torque (N-m)

Visc. 7 Pa.s

Visc 15 Pa.s

Visc. 27.25 Pa.s

Figure 4.10: Time evolution of torque in anchor mixer for various viscosities.

In this situation, a sufficiency of the species transport equation can be spoiled by the

effect of so-called numerical diffusion, i.e. the artificially generated smoothing of a

tracer profile due to errors from the discretisation, and the numerical method can be

much stronger than the true physical diffusion. In fact, almost all of the mixing that

is observed in a numerical simulation may then be artificial. To show the intensity of

numerical diffusion we performed simple test case in a lid driven cavity flow. Fig. 4.11

shows the influence of solver settings on numerical diffusion i.e., the so called dis-

cretization. The left hand side figure shows the tracer profile under the approximation

of 1st order discretization (UDS) and yields produces numerical diffusion due to Taylor

series truncation error [Eq. 3.4]. The figure on the right hand side depicts the tracer

profile under the approximation of 2nd order discretization (CDS) and numerically

less diffusive because the approximation is second order [Eq. 3.5]. One way to reduce

this numerical diffusion is to refine the mesh as fine as possible as shown in Fig. 4.12.

But the problem with mesh refinement is that it is numerically very expensive and the

numerical diffusion is only minimized and not keep off.

A way to avoid this problem is to replace the continuous tracer concentration by

60 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

Figure 4.11: Influence of solver settings on numerical diffusion. 1st order discretization

(left) and 2nd order discretization (right).

Figure 4.12: Influence of computational mesh refinement on numerical diffusion.

Coarse mesh (left) and fine mesh (right).

4.3. PARTICLE TRACKING 61

1 2 3

Figure 4.13: Contour plot of number concentrations (right) obtained from the number

of tracked particles inside compartments.

a number concentration obtained from Lagrangian (i.e. inertia free) particles that are

tracked during the simulation. This approach does not suffer from artificial diffusion,

since the position of tracer particles can be resolved with sub-grid-scale accuracy and

the velocity field at these particle positions can be obtained by interpolation from the

grid values. The artistic view of contour plot as a result of number concentrations ob-

tained from the number of tracked particles inside compartments is shown in Fig. 4.13.

Within this approach, initial positions x1,...,xnfor a certain number of nparticles

are chosen-typically with random positions inside an appropriate subregion such that

they form a kind of blob. Then the particles trajectories are computed as the solution

of the dynamical system:

x(t) = v(t, x(t)) ,x(0) = xi.(4.22)

In addition, a certain number mof compartments is defined, where the compartments

form a subdivision of the mixing region, i.e.

[

k=1

Vk,with disjoint Vk.(4.23)

Let Xk(t)denote the number of particles that are inside Vkat time t. Then the

62 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

5 30

0.05

0.1

0.15

0.2

0.25

0.3

Figure 4.14: Decay of normalised variance under application of the Baker map. Effect

of numbers of particles (P) and compartments (C): 100P, 64C (red); 800P, 4C (blue);

3200P, 64C (black).

normalized number concentrations ck(t)are defined as:

ck(t) = Xk(t)

|Vk|/n

|V|,(4.24)

which reduces to:

ck(t) = m

nXk(t),(4.25)

in the case when all compartments Vkhave the same volume. The basic question

within this approach is how many particles and how many compartments have to be

used to obtain a reliable mixing time from the decay of variance of the normalized

number concentrations? A particular time evolution of the variance is of course only

a single realisation of a stochastic process. Hence any quantity obtained from it is a

random variable that can only be characterized by means of its statistical properties.

Furthermore, a state of complete homogenization has a defined meaning only in the

statistical sense. Especially, such a state will show a strictly positive variance. There-

fore, in order to be able to compute a reasonable mixing time, this "‘residual variance"’

4.3. PARTICLE TRACKING 63

Stretching FoldingStretching StretchingFolding

Folding

Stretching FoldingStretching StretchingFolding

Folding

Figure 4.15: Illustration of Baker map

Figure 4.16: Typical particles trajectories in Lid driven cavity.

σ2

Rhas to be much smaller than the threshold s2(0) from above. This is illustarted

in Fig. 4.14, which shows the decay of variance under successive applications of the

Baker map to a set of nparticle positions that are initially close together. Let us note

in passing that the Baker map was introduced in the 1930s in [66], probably inspired

by Birkhoff. Nowadays, the Baker map is a well-known prototype map that models a

simple mixing process composed of stretching and folding. Its definition is implicitly

given in Fig. 4.15. In Fig. 4.14, red dots correspond to the variance obtained from

100 particles and 64 compartments. Evidently, the residual variance is too large to

reach the threshold for = 0.05. In case of 800 particles, but only 4 compartments,

the residual variance is below the critical value, but the stochastic deviations of these

variances from their mean, i.e. the variance of variances, is too large. Reliable results

are obtained, e.g., in case of 3200 particles and 64 compartments as shown by the

black dots.

64 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

Figure 4.17: Schematic representation of anchor mixer with hybrid (combination of

tetra and hexa) mesh. The mesh contains about 100000 control volumes and is gen-

erated using Star-Design.

Figure 4.18: Lagrangian particles trajectories in anchor mixer.

4.3. PARTICLE TRACKING 65

11 10

160

127

147 16.5

Figure 4.19: Schematic representation of anchor mixer with dimensional configurations

(in mm) and it is divided into 4 compartments vertically for particle tracking. A blob

of particles also can be seen in the figure.

Figure 4.20: Schematic representation of anchor mixer with division of 16 compart-

ments in radial and circumferential cross section.

66 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

A similar problem concerning residual variances occurs if mixing times are deter-

mined from experimental measurements. In this case σ2

Rhas been calculated under the

assumption that the tracer particles are evenly distributed and that the resulting con-

centrations in individual probes are stochastically independent. The latter assumption

is not reasonable if the complete volume is samples by means of all compartments,

since then:

k=1

pkck(t) = 1 with pk=|Vk|

|V|,(4.26)

i.e. the random variables ck(at a fixed time t) are not independent. In fact, if the

particles are evenly distributed in Vand if Xkdenotes the number of particles inside

Vk, then their joint distribution is given by the multinomial distribution:

P(X1=l1, ....., Xm=lm) = n!

l1!...lm!pl1

1....plm

m,where n=

k=1

lk.(4.27)

We restrict our attention to compartments of equal volume, in which case the above

equation simplifies to:

P(X1=l1, ....., Xm=lm) = n!

l1!...lm!

mn,where n=

k=1

lk.(4.28)

Then, by an elementary calculation yields:

Var[Xk] = n1

m(1 −1

m),hence σ2

R=Var[ck] = m

n(1 −1

m).(4.29)

If, initially, the particles are evenly distributed in lof the mcompartments, the initial

variance is:

σ2

0(X) = n2

m(1 −l

m)1

l.(4.30)

This leads to a relative residual variance of:

σ2

0(c)=Var[Xk]

σ2

0(X)=n

m(1 −1

m)1

m(1 −l

m)1

=(m−1)l

n(m−l)≈l

nfor lm . (4.31)

The estimate Eq. 4.31 can now be used to calculate the number of particles needed to

compute a reliable mixing time t1−

Mfor a given level  > 0. Indeed, this number nhas

to chosen such that:

γ=l

n 1,(4.32)

4.4. RESULTS AND DISCUSSIONS 67

say γ=0.1. As an example, for =0.05 and l=1, at least 200 particles should be

used, but 2002= 40000 particles are required if the mixing time is computed from

standard deviations instead of variances. Within the computations, the true variance

is estimated by means of:

s2(X) = 1

k=1

(Xk−n

m)2,resp. s2(c) = 1

k=1

(ck−1)2.(4.33)

To avoid wrong mixing times caused by too much noise contained in the computed

variance, one should also check for:

Var[s2(c)]

σ2

≤γ , (4.34)

say. Again by elementary calculations it follows that:

Var[s2(c)]

σ2

=Var[s2(X)]

E[s2(X)] =Var[s2(X)]

Var[X]=2

m−1(1 −1

n)≈2

m−1(4.35)

which determines the number of compartments needed. Surprisingly, quite few com-

partments are sufficient to keep the noise at a small level.

4.4 Results and Discussions

The above developed method has been evaluated using the numerical investigation of

one of the oldest and yet widely used mixer configuration which is shown in Fig. 4.19.

It consists of a flat bottomed vessel of 147 mm height and a diameter of 160 mm.

The tank is equipped with anchor impeller of width 127 mm and height equal to the

vessel height of 140 mm. Other important dimensions such as clearance between the

wall and anchor, blade width and thickness are shown in Fig. 4.19. Unstructured

hybrid (combined hexa and tetra type) mesh consists of approximately 70,000 control

volumes was generated using Star-Design software. The working liquid is a Newtonian

fluid with material properties density: ρ=1000 kg/m3 and viscosity: µ= 5 Pa s. A

speed of 50 rotations per minute is set for the rotating reference frame. Since the

boundaries of the rotating domain is assumed to be rotating witht the same angular

speed, an opposite spin of the same angular speed is defined on the inner walls to make

68 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

them stationary. Reynolds number of 2.7 is held which is correspond to laminar flow

condition. The time dependent transport equations [Eq. 3.12] are solved in Cartesian

co-ordinates for laminar flow. The influence of grid dependency is investigated at 3

different grid sizes and the observation is that the solution has a very little effect

based on the velocity field comparison. A compramise in grid size is arrived based

on accuracy of the calculation and numerical efforts needed. A second order accurate

spatial discretization [Eq. 3.5] scheme is used to solve convective terms. The following

section discusses the calculation of intensity of segregation.

4.4.1 Intensity of Segregation

In case of transient species distribution c(t, x), the standard deviation σ=σ(t)and,

hence, the variation coefficient will also depend on time, although the mean value stays

constant if no flux of tracer across the domain boundaries occurs. The evolution of σ(t)

indicates how the quality of mixing changes in the course of time occurs. The evolution

of the homogeneity of a scalar field at a given time, it cannot be directly employed for

the definition of a mixing time or a mixing length. This is due to the fact that it is not

normalized to the range 0...1, say, but can attain any value between zero and infinity.

Hence, while an ideally homogeneous mixture has a segregation coefficient of zero, it is

not clear which threshold should be reached during a mixing process. An appropriate

normalization can be derived more easily in case of non-reactive mixing and for an

isolated mixing region. As mentioned above, the expectation µstays constant then,

hence only the variance σ2needs to be normalized. This leads to:

I=σ2

σ2

max

,(4.36)

where σ2

max denotes the maximum variance. These considerations result in Danckw-

ert’s intensity of segregation [16], defined by:

IS=σ2

µ(1 −µ).(4.37)

Based on Danckwert’s intensity of segregation, measures for the intensity of mixing

can be derived and a large number of related measures which are similar in the sense

4.4. RESULTS AND DISCUSSIONS 69

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 10

Time (s)

Normalised Variance

Figure 4.21: Time evolution of normalised variance for anchor mixer. Effect of numbers

of particles (P) while compartments (C) kept constant: 100P, 64C (red); 1600P, 64C

(blue).

that they depend on statistical parameters are reviewed in [12]. A typical example,

which is recommended there, is:

IM= (1 −qIS) = 1 −σ

σmax

.(4.38)

Since ISis normalized, this measure is zero for completely segregated mixtures and

attains the value one in the homogeneously mixed case.

In order to assess the quality of mixing in stirred vessels, the significant task was

to compute the tracer particles trajectories and obtain the data for variance evolution.

To achieve this, a mass-less Lagrangian particles of 1600 are randomly placed initially

in one compartment (sub-volume) such that they form a blob of particles as shown in

Fig. 4.19 and Fig. 4.20. This approach does not suffer from artificial diffusion, since

the position of tracer particles can be resolved with the sub-grid scale accuracy and the

velocity field at these particle positions can be obtained by interpolation from the grid

values as mentioned earlier. Notice that the co-ordinates (position) of the particles

70 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

0.00

0.04

0.08

0.12

0.16

0.20

60 70 80 90 100

Time (s)

Normalised varianc

Figure 4.22: Time evolution of normalised variance for anchor mixer. Effect of numbers

of particles (P) while compartments (C) kept constant: 100P, 64C (red); 1600P, 64C

(blue).

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 10

Time (s)

Normalised varianc

Figure 4.23: Time evolution of normalised variance for anchor mixer. 1600 Particles

are initially placed in one compartment (blue), multiple compartments (red).

4.4. RESULTS AND DISCUSSIONS 71

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 10

Time (s)

Normalised Variance

Figure 4.24: Time evolution of normalised variance for anchor mixer. Effect of numbers

of compartments (C) while particles (P) kept constant: 1600P, 4C (red); 1600P, 64C

(blue).

are in rotating reference frame and it is necessary to transform them to laboratory

frame. The particle number information is recorded during the simulation process

with the aid of user defined subroutine written in FORTRAN 77 and provided to Star-

CD. During the process of particle tracking, the number concentration is calculated

according to Eq. 4.24. The evolution of variance can then be computed based on

number concentration information as defined in Eq. 4.33.

Time evolution of normalised variance (is also known as intensity of segregation,

obtained by dividing computed variance with the initial variance) can be seen in

Fig. 4.21. The value on Y-axis varies between one and zero. Intensity of segregation

(Is) equal to 1 means complete segregation and a value of 0 stands for complete

homogenisation. To achieve complete homogenisation (Is=0) in numerical studies

is very difficult even with the present day computers. In Fig. 4.21, the blue curve

correspond to Iswith 1600 particles and the red curve correspond to Iswith 100

particles. In both the cases, particles are initially placed only in one compartment. A

72 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

sharp decrease is observed at the beginning of the calculation in these curves which

is attributed to the domination of the flow due to kinematical mixing. When time

proceeds, a condition of unvarying in the evolution curves occur which is due to the

fact that almost even distribution of the particles in the whole domain is obtained.

The one of the important observations from this figure is that the red curve possesses

more noise than the blue curve (Fig. 4.22), owing to insufficient number of particles.

In the case of anchor mixer, it is well known fact that the velocity gradients are poor

in vertical direction. Hence 1% of initial variance is achieved only after expensive

numerical simulation.

In Fig. 4.23 we illustrate the effect of number of compartments where particles are

initially placed. It shows two curves with different colours, the one in blue is the result

of time evolution of Iswhen the particles are initially placed in one compartment and

the other in red colour is the result of Iswhen the particles are placed initially in

multiple compartments. Remarkable difference can be seen in the magnitude of Is

between these two curves. If particles are initially placed in more than one compart-

ment, there is a possibility of attaining Isgreater than one, since after certain time

elapses particles may be advected to the compartments which are less than the com-

partments in which particles are placed initially. The other important observation in

the current investigations is that the number of compartments required for reasonable

mixing quality. This phenomenon is depicted in Fig. 4.24. If the total computational

domain is divided into few compartments, the residual variance is too large and also

the fluctuations in Isare too high even with more particles.

4.4.2 Mixing Time

A general definition of σ2

max is problematic, since a universal maximum value for the

molar concentration does not exist. One possible adaptation to molar concentrations

or, more generally, intensities is to let σ2

max be the variance of a totally segregated

distribution composed of the same amount of substance or the same total intensity

and the same maximum value cmax. Again, this variance is independent of the concrete

4.4. RESULTS AND DISCUSSIONS 73

structure of the segregated field and leads to σ2

max =µ(cmax −µ). This reference value

allows for a normalization of the variance of the distribution based on the scalar field

itself., i.e. without need for an additional external value. This is useful for assessing

a species distribution by itself, but, since cmax will usually change with time or space

it cannot be employed for defining a mixing time or a mixing length. On the other

hand, if a mixing time or length is to be defined, a mixing process instead of a single

mixture is to be evaluated. In this case, mixing evolves in the course of time or along

a spatial direction and, hence, the variance σ2

0at the initial time or the entrance of

the mixing channel can be used instead of σ2

max. For instance, a mixing time tmix is

obtained as the first time t0such that:

s2(t)≤s2(0) for all t≥t0,(4.39)

where (with 0<<1) is a given fraction; typical values are =0.05 or =0.01. This

criterion corresponds to the condition σ2/σ2

0≤or to IS≤with ISfrom Eq. 4.37 in

case of a completely segregated initial state. Of course, these mixing times strongly

depend on the choice of . As a variant of the later definition above, a mixing time is

also defined as the first time from which on IMfrom Eq. 4.38 stays above 1-.

The determination of mixing times from numerical tracer experiments with the

help of CFD simulations can be difficult, since a value of IM=0.99, say, will often be

only reached after long and expensive calculations. This is especially true for highly

resolved simulations which need to be used in case of high Schmidt numbers. On

the other hand, there are several references, apparently starting with [40] in which an

exponential asymptotic behavior of the type:

σ2(t)∼exp(−t/τ∞

mix)for large t , (4.40)

with an assymptotic mixing time τ∞

mix is assumed. In case of a continuous mixing

process, this type of behaviour with more time replaces by axial position is employed

in [31] to define a mixing length. While Eq. 4.40 is often used without justification,

there is increasing mathematical evidence for Eq. 4.40 to be valid in several situations

[44] and [75]. In situations where Eq. 4.40 holds, this can of course be used to reduce

74 CHAPTER 4. CFD-ANALYSIS OF ANCHOR MIXER

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 10

Time (s)

Normalised Varaince

Figure 4.25: Time evolution of normalised variance for anchor mixer with 1600 parti-

cles and 64 compartments.

the computational effort needed for numerical mixing time calculation. Indeed, it then

suffices to run a numerical tracer experiment until a linear dependence of logσversus

time is reached and the mixing time can then be obtained by extrapolation [29].

In the case anchor mixer, the mixing time is computed based on the procedure

mentioned above and let us consider the value of = 0.05. As it can be seen from

Fig. 4.25, normalised intensity never reaches a value of = 0.05. If we consider

= 0.1, normalised intensity achieves certain minimum but again increases. Hence the

appropriate definition of mixing time in anchor mixer is difficult and even after a long

and expensive computations. These observations reveal the poor mixing characteristics

of anchor mixer. It is well known that anchor mixer is good for heat transfer between

wall and the fluid rather than its mixing capabilities. Numerical investigation of

mixing in a specific kneader element is presented in the next chapter.

Chapter 5

Mixing in a Kneader Element

5.1 Introduction

Synthesis of highly viscous materials such as polymers, rubber and food etc requires

suitable machines and screw extruders are well known for this purpose. However, screw

extruders are no more efficient when longer residence time and little axial dispersion is

needed. Hence, several investigations were made in the early 20th century to modify

screw extruders to achieve longer residence time and little axial dispersion. One of

the such maschines is kneader reactor which is used for mixing of highly viscous

liquids especially paste like materials and these maschine has a construction similar

to screw extruders. During this time the idea of using extrudes with pins on a barrel,

reciprocating screws and self-cleaning action evolved. The purpose of pins [3, 84]

on barrel is to agitate and knead the material being processed to remove air from

it. The reciprocating action of screws yield better injection molding for highly viscous

polymers. The term self cleaning means here every surface which is in contact with the

reacting product is mechanically wiped by another surface with a well defined closed

distance. The principle of rotating and oscillating action of single shaft Kokneader

is shown in Fig. 5.1. Keeping the screws clean to avoid collection of dirt and sticky

nature of polymers to the kneading pins and barrel can be achieved by self cleaning

mechanism [Fig. 5.2] of the extruder.

76 CHAPTER 5. MIXING IN A KNEADER ELEMENT

Figure 5.1: Single shaft kokneader with oscillating and rotating mechanism; 1. Barrel,

2. Kneading pin, 3. Kneading element, 4. Shaft, 5. Oscillation and 6. Rotation.

Heinz List [43] invented kneader by combining the better kneading, injection

and self cleaning mechanisms. Further developments were intensely took place be-

tween 1960-70 to further modify kneaders for various applications. For more detailed

overview on the literature related to kneader development can be found in [83] and the

references given there. Few experimental investigations on buss kneader characteris-

tics were done between 1955-95. However, the first basic and detailed experimental

study is reported by [20]. They used silicone oil (µ= 1 Pas) and a parafanic oil (0.2

Pas) under isothermal conditions at room temperature to measure throughput versus

screw speed, pressure and fill lengths. Little later, Lyu and White [45, 46] published

their findings on fill factor and temperature profiles, melting conditions and residence

time distribution as a function of screw configuration and processing conditions for

List/Buss kneader. More recent studies [76,77] concerning Buss kneader as a polymer-

ization reactor shows the good mixing properties of kneader in radial as well as axial

directions for exothermal reactions. Simulations of fluid flow in a Buss kneader started

in the late 80s to model flow and pressure field distributions. A detailed investigations

on total machine characteristics of Buss kneader were published by Elemans and Mei-

jer [20], Lyu and White [45–48]. A 3-D computational study reported by Mehranpour

et al. [53] to predict the velocity field in the conveying element of a Ko-kneader shows

5.1. INTRODUCTION 77

Figure 5.2: Schematic of the self-cleaning mechanism of kneader by movements of the

pins [83].

that the reciprocating action of Kokneader enhances the mixing performance by means

of periodically changing the flow field and shear rate distribution.

In addition to single shaft, kneaders with double shaft are also gained much at-

tention in the recent years. The advantage being large length to diameter ratio and

shows clear analogous to twin screw extruder. A typical configuration of double shaft

kneader is shown in Fig. 5.3 which has large length to diameter ratio and also the

incident angle of the two shafts. The angle formed by the surface of each kneading

element with the shaft of the rotor towards the bottom of the tank is designated

by αwhere as the angle formed by the surface of each kneading element with the

shaft of the rotor against the bottom of the tank is designated as β. Similar to twin

screw extruder, the flexibility of double shaft kneader is also that it can be corotat-

ing or counter rotating. The experimental works are mainly limited to determination

of integral, characteristic features of the whole reactor, such as mean residence time

or mean degree of filling [5, 63]. The experimental investigation of residence time in

kneader, single as well as twin screw extruders are carried out by several researchers

78 CHAPTER 5. MIXING IN A KNEADER ELEMENT

using different tracer materials, for example using aluminium flakes [49,67], coloured

materials [5,15,62,69,73], iron powder [63], salts [17,19,78,79] and light emission of

flourescent tracer [15]. However, local phenomena and fundamental mechanisms are

not studied. The recent master thesis [65] reports on experimental investigation of

double shaft intermeshing kneader. Their results include investigation of power dis-

sipation, fill level, residence time distribution (RTD) and heat transfer for the better

understanding of the basic mechanisms in the kneader.

For the description of the process, simple cumulative models such as the cascade

model [48,54], dimensionless numbers relations or even regression approaches [51] are

employed but cannot contribute to a thorough understanding. Only very little work

is based on continuum mechanical models. These are mainly devoted to twin screw

extruders. In that case, flow geometry make the use of lubrication theory relevant to

fully filled parts, resulting in 1D-models for RTD, average temperature and pressure.

For simplified materials rheology, these models, together with specific optimization

algorithms, can help to optimize compounding processes up to a limited level only

[25, 42]. Recent progresses using 3D local approaches have also been made [7] using

fluid-structure interaction techniques [28]. Admittedly, according to the knowledge of

the author of the present thesis, there are no numerical investigations reported in the

literature on the investigation of double shaft kneaders. In particular, essential aspects

of kneaders such as interpreting geometries of the rotating shafts, the appearance of

free surfaces and, above all, an adequate modelling of mixing and mass transfer are still

missing. The reason is that the limited capability of numerical simulation techniques

to handle intermeshing or ovelapping grid problems. Basic numerical investigation of

a particular unit cell of a kneader element is presented in this chapter. The geometry

used in the current investigation is a cylindrical stirred tank equipped with complicated

kneading disks and pins as shown in Fig. 5.4. This geometrical configuration contains

kneading disks attached to the shaft and kneading pins mounted on the kneading discs

for the purpose of cleaning the barrel if any material is sticked to it especially during

highly viscous polymer mixing. Simulation procedure in this stirred tank is explained

5.2. SIMULATION 79

Figure 5.3: Two shaft intermeshing Reacom 60L kneader reactor.

in the following section.

5.2 Simulation

For the computations presented here, a general purpose finite-volume based commer-

cial CFD package from CD-Adapco has been used. This package includes tools for

parametric geometry definition, automatic mesh generation, flow-field solution and

postprocessing the results.

5.2.1 Flow Computation

Time dependant simulations were performed for the flow created by kneading element

in a tank. The tank itself is a vertical, cylindrical vessel with a diameter of 130 mm

and height of 90 mm. The concentrically located shaft extends the entire height of the

tank and has a diameter of 60 mm. Each of three kneading disks are spaced at 120°

in circumferential direction and each of three kneading disks in vertical direction are

spaced at 40° each in vertical direction corresponding to equally sized nine kneading

disks. Kneading pins are appropriately mounted on the kneading disks without losing

80 CHAPTER 5. MIXING IN A KNEADER ELEMENT

symmetry. It should be noted that the computational model of this system represents

the exact geometry including the small and angled pins. The unstructured mesh is

generated using Star-Design and contains approximately 100,000 hybrid (hexa and

tetra) computational cells. The mesh is exported to Pro-Star where the parameters of

simulation are provided. The properties of the working liquid in this case are constant

viscosity (µ) of 1 Pas and constant density (ρ) of 1000 kg/m3. The laminar, incom-

pressible and time dependant momentum conservation equation in rotating reference

frame [Eq. 3.12]:

∂u

∂t + (u· ∇)u+Ω×(Ω×r) + 2Ω×u=−1

ρ∇p+ν∇2u+g,(5.1)

is solved, where gis the acceleration due to gravity. The continuity equaion (∇·u= 0)

is invariant regardless of the reference frame. The whole computational domain is

rotated (solid-body rotation) with an angular velocity of 60 rotations per minute.

The velocity boundary conditions for all cells on the kneading disks, pins together with

shaft are set to zero by defining equal and opposite magnitude of angular velocity used

for solid-body rotation. The velocity boundary conditions for all computational cells

on the cylindrical vessel surface, top and bottom walls are set to correspond to the

solid-body rotation. The initial velocities of the fluid cells are assumed to correspond

to the same solid-body rotation as the tank walls. A time step of 0.0005 sec is set and

several time steps were carried out resulting adequate rotations of the impeller.

5.2.2 Particle Tracking

In the present study, mass-less fluid particles are introduced in the flow is followed

using a Lagrangian particle tracking method in order to quantify the mixing quality.

A randomly distributed 1600 particles are placed in small volume in the middle of

the vessel at the beginning of the simulation process. As discussed earlier, particle

tracking approach avoids the introduction of numerical diffusion which results if a

scalar variable is tracked, which confuses the actual mixing behaviour of equipment.

The movement of the particle tracers in the flow is tracked by the integration of the

5.3. RESULTS AND DISCUSSIONS 81

(a) (b)

Figure 5.4: Schematic of (a) unit cell of a specific kneader and (b) the top view of

kneader element (CAD data from BASF AG, Ludwigshafen).

vector equation of motion of each particle [Eq. 4.22]. A rebound condition given to

the particles to avoid trajectories being trapped near the walls where the local velocity

is close to zero.

5.3 Results and Discussions

Simulations are performed in a kneader element as explained in the previous sections

and the post processing of the results are discussed below.

5.3.1 Flow Patterns

In order to understand the flow inside the kneading element, velocity magnitude sec-

tion plots are discussed here. Fig. 5.5(a) shows the velocity contour plot in horizontal

cross section. Red colored regions denote high velocity magnitudes, blue colored re-

gions denote low velocities while the regions between red and blue color denotes the

intermediate velocity magnitudes. The highest velocities are observed near the tip of

the kneading disk, while velocities are lowest near the walls. Fig. 5.5(b) shows the

82 CHAPTER 5. MIXING IN A KNEADER ELEMENT

velocity vectors in 3-D horizontal cross section made at the middle of the vessel. The

vectors point in the direction of liquid velocity at the point where they originate. The

length of the vectors is proportional to the magnitude of the liquid velocity. The flow

is weak near the walls and strong around the kneading disks. It is also found in both

the plots that the regions with low velocities are minor and hence the high mixing

behaviour of kneading element is justified.

(a) (b)

Figure 5.5: Horizontal section plots of velocity fields. (a) Contour plot and (b) vector

plot.

5.3.2 Mixing Quality

It follows the same method as discussed in previous chapter. In order to quantify

the homogeneity of a mixture, a statistical analysis of the concentration in samples

from the mixture, which is based on Danckwert’s intensity of segregation concept [16]

can be used. The intensity of segregation approach is based on the variance of the

concentration at different regions in space with respect to the mean concentration. In

the case of massless tracer particles, a number based variance is calculated by dividing

whole computational domain into small compartments. The intensity of segregation

is then variance divided by its maximum value.

An intensity value of 1 corresponds to complete segregation and a value of 0 in-

5.3. RESULTS AND DISCUSSIONS 83

0.0

0.2

0.4

0.6

0.8

1.0

0 3 6 9 12

Time (s)

Normalised Variance

Figure 5.6: Time evolution of normalised variance for kneader element with 1600

particles and 64 compartments.

dicates complete mixing. To determine the mixing quality, the most important task

was to compute the tracer particles trajectories. This has been usually done based

on the velocity integration over time. To do this, total computational domain was

divided into 64 compartments (sub-volumes) of approximately equal volume and 1600

particles are randomly distributed initially in one of the compartments. Then the

variance and also the intensity of segregation is computed. The results obtained are

presented in Fig. 5.6 where the evolution of intensity of segregation is plotted against

time. The value of Isstarts at 1 where the particles are completely segregated and

reaches close to zero corresponding to the optimal dispersion. It can be seen that the

intensity of segregation is decreasing in unordered manner with random periods at the

beginning of the simulation. This is due to the fact that the tracer particles are still

need to undergo significant stretching. It can also be observed that the level off of

the intensity of segregation curve occurs at time 8s. This time can be defined as the

mixing time as it indicates when the tracer concentration stabilizes. The intensity of

segregation values after reaching this mixing time is = 0.05. However, we still can

go below this value to about = 0.007, but again the intensity of segregation starts

84 CHAPTER 5. MIXING IN A KNEADER ELEMENT

to increase.

5.4 Mapping Method

In the following section, mapping matrix approach is described for the simulation of

mixing in a specific kneader element. The method is based on a spatial discretization

of the locally averaged concentration of fluid components in the mixture. A distribu-

tion matrix, that describes the changes in a component concentration, is composed.

This mapping matrix approach makes it possible to rapidly predict the mixing per-

formances. Introduction and basic principle of this approach is given in the following

subsection.

5.4.1 Introduction

As already mentioned in the previous chapters, mixing of fluids is a topic of signif-

icant interest, because of wide applications in process industries. One of the simple

but important class of mixing phenomena is laminar mixing. Inspite of considerable

achievements in the understanding of its mechanisms, numerical mixing remains com-

putationally expensive and requires special methods. Further, determination of the

accurate velocity field is a very important step. Once the velocity field is known, a

variety of techniques is available to study mixing based on the tracking of deforming

individual fluid volumes such as front capturing and front tracking. Front capturing

technique describe the advection of a special material function and accordingly use

post-processing techniques to restore the interface shape. In contrast to this, front

tracking techniques rely on an explicit description of the interface shape for which

an auxiliary surface mesh is used. This approach yields an accurate description of

strongly deforming fluid volumes and is capable to deal with the mixing of different

fluids (distributive mixing). But the amount of data produced increases exponentially

in case of efficient mixing flows. Most spatially bound mixing flows exhibit temporal

or spatial periodicity, dynamical system methods [58] such as Poincaré maps, peri-

5.4. MAPPING METHOD 85

odic points and analysis of their manifolds are useful. Poincaré maps help to reveal

zones of chaotic mixing and regular motion where as periodic points are regarded as

chaos or regularity. However, they do not provide information on the rate of mixing

and concentration distribution. Numerical tracer experiments are of increasing impor-

tance to analyze mixing, but the accurate solution is spoiled by numerical diffusion.

Lagrangian particle tracking can give the particle distribution concentration during

the course of the mixing, but requires lot of computational memory and time. A way

to avoid these problems is using of mapping matrix method.

Since chaotic mixing of viscous liquids in laminar flows (ex. kneaders, static and

micro mixers) is usually based on the situation where the Baker’s transformation is

applied a number of times on a specified volume of material, distribution of material

in such flows can be handled quite well by distribution matrix method. The mapping

method is essentially the transport of a conservative quantity by means of mapping

matrix, describing the transport of fluid from an initial cross section to final one in

the case of space periodic flows or from an initial time to final time in the case of time

periodic flow (in the present case). The advantage of using mapping method is that

it requires just one-time computation of the deformation induced by the flow during

fixed flow in time ∆tand the effect of flow for any number of combination of cycles

can then be evaluated by a repeated multiplication of the distribution matrix with a

prescribed initial concentration distribution vector. Since these multiplications take

only a few CPU seconds, this brings a huge benefit over conventional tracking method

in which the tracking is repeated from the first to last period to analyze mixing which

results in lot of computational time. The basics and the mechanism of the mapping

method is described below.

5.4.2 Basic formulation

Laminar, highly viscous liquid mixing is typically based on a repetetion of the same

motion characterized by three simple steps: streching, cutting and folding a specified

volume of material. The dustribution of the material in such flows can be handled

86 CHAPTER 5. MIXING IN A KNEADER ELEMENT

AtD

Figure 5.7: Schematic of the mapping matrix method: stretching of fluid element.

quite well by the use of mapping matrix methods as originally proposed by [71] and

was redefined in a computationally efficient way by [37]. In this section we describe the

formulation of mapping method by considering two cases. One is based on the tracking

of boundaries of a sub-volumes which cover the entire fluid volume and the other is

based on the tracking of particles in the compartments of the complete computational

domain.

Case-1: Let us consider an arbitrary fluid region Vis divided into a number of

sub-regions Vkthat covers the entire fluid volume. The boundaries of the sub-voulmes

are tracked from t=t0to t=t0+ ∆t. The two grids are superimposed (Fig. 5.7) and

the intersection Aij is computed as:

Aij =Z

Vj(t=t0+∆t)TVi(t=t0)

dS/ Z

Vj(t=t0)

dS , (5.2)

where Aij is the fraction of sub-volume Vjat t=t0, that is tracked to t=t0+ ∆tand

found in the sub-volume Vi. These intersections are stored in a matrix called mapping

matrix or distribution matrix A. To obtain the intersection points of the distribution

matrix, the total cross section of the flow domain is sub-divided into a large number of

discrete cells (N)of identical cells. During the flow, the material from a donor cell is

advected to different recipient cells. The fraction of material that is transferred from

the donor cell to a recipient cell gives the distribution coefficient of the donor cell with

5.4. MAPPING METHOD 87

t'j

Figure 5.8: Schematic of the mapping matrix method: particle distribution.

respect to the recipient cell. Therefore, in total Ncells from mapping matrix Aof

the order NXN. For further information on this method please refere to [37]. But

tracking all interfaces of all Ncells during a flow over a time ∆tis cumbersome to

track interfaces experiencing complicated deformation patterns. Thus, an alternative

approach that is much simpler is presented below.

Case-2: Similar manner is followed as above to generate a mapping matrix using

particles number concentrations to compute variance. A detailed formulation of the

mapping method based on the particle number concentration is given in [68]. In this

case, the coefficients of the mapping matrix Aare calculated in the formulation given

here. A schematic representation of the particle distribution in the sub-volumes are

shown in Fig. 5.8. Particles inside all sub-volumes are tracked to approximate the

coefficients of A. The particles are uniformly distributed in the sub-volumes. If the

number of particles in sub-volume Vjis Bjat t=t0and the number of particles found

after tracking in the sub-volume Viis Bij at t=t0+ ∆t, then the mapping coefficient

Aij is computed as:

Aij =Bij

.(5.3)

In general, the coefficient Aij is the measure of the fraction of total flux of the sub-

volume Vjdonated to the sub-volume Vi. Since the majority of the efficient mixers

show periodicity either in space or time, the same matrix is repeated with itself to

obtain Anin order to account for the periodicity of the flow as suggested by [71].

88 CHAPTER 5. MIXING IN A KNEADER ELEMENT

By studying the properties of An, mixing can be analyzed by its evolution in time

or space. However, this procedure is not so attractive and feasible especially in 3D-

situation where An(dense matrix) generates large data and hence require huge storage

capacity. Instead of analyzing the evolution of the matrix, the concentration Cof a

fluid in each sub-volume is computed as proposed by [24, 37,38]. Thus, alternatively

Cnis computed by the sequence for i=1 to n:

Ci=ACi−1.(5.4)

This procedure eliminates the necessity to compute Anand therefore, much cheaper in

numbers of operations as well as computer memory. The mapping matrix calculations

are easily parallelized [23].

5.4.3 Application to a Kneader Element

Here, we consider unit cell of a single shaft specific kneader element for the application

of mapping method. First, computations are performed using 1600 mass-less particles

placed initially in one of the 64 sub-volumes (compartments) and collected information

regarding the particle distribution in each compartment at regular time intervals. Once

all the compartments are occupied with approximately equal number of particles,

simulations restarted and tracked the trajectories of each particle. Each particle can

be identified with particle ID number. This means, for example a particle with number

1 in compartment 1 at t=t0, advected to some other compartment after t=t0+ ∆t.

Here, ∆tcan be the time required for one or and one-third rotation of the impeller.

Likewise, positions of all 1600 particles are identified after one rotation and stored

them in a matrix of an order 64X64. The stored data contains one starting position

and the end position for each numbered particle. Matrix coefficients are created by

comparison of starting and ending position. For every particle starting in compartment

mand reaching the compartment n, the coefficient in the mapping matrix is raised by

one. Then, each row of the matrix is normalized with the total number of particles

in the complete row. Multiplication of the resultant matrix with the matrix vector

5.4. MAPPING METHOD 89

1.0E-04

1.1E-04

1.2E-04

1.3E-04

1.4E-04

1.5E-04

1.6E-04

6 9 12 15 18Time (s)

Variance

Mapping

Tracking

Figure 5.9: Time evolution of variance in a kneader element: Comparison of simulated

data with data obtained from mapping matrix.

containing single distribution (at t= 0, single column matrix) yields mapping matrix.

First case, time evolution of variance is computed with the use of the resulting mapping

matrix at this particular time step. This procedure is repeated without using CFD

software further to obtain the variance at each time step for up to 10 rotations. On

the other hand, CFD simulation is performed for the same case file and computed

the time evolution of variance for each rotation of the impeller. Fig. 5.9 verifies the

variance values computed using mapping matrix method and simulation using CFD

software. The solid line represent the variance values obtained using tracking where

as the points correspond to mapping. In this case, mapping matrix at corresponding

rotation is multiplied with the matrix from the previous rotation.

Next, we consider the mapping matrix obtained only after one rotation and used

to compute the particle distribution with base matrix with a vector containing the

initial distribution. The distribution after nrotations can be obtained by multiplying

the initial distribution vector with the matrix for ntimes. Although one can use

every desired initial position of the particles, simulation data by tracking is needed for

comparison. The comparison of the time evolution of normalised variance computed

90 CHAPTER 5. MIXING IN A KNEADER ELEMENT

0.2

0.4

0.6

0.8

1.2

0 2 4 6 8 10 12 14 1

Time (s)

Normalised Variance

Tracking

Mapping

Figure 5.10: Time evolution of normalised variance in a kneader element. Comparison

of simulated data with data obtained from mapping matrix after every rotation of the

impeller.

0.2

0.4

0.6

0.8

1.2

0 2 4 6 8 10 12 14 1

Time (s)

Normalised Variance

Mapping

Tracking

Figure 5.11: Time evolution of normalised variance in a kneader element. Compari-

son of simulated data with data obtained from mapping matrix after every one-third

rotation of the impeller.

5.4. MAPPING METHOD 91

based on mapping method and CFD simulation is shown in Fig. 5.10. Each marked

point in this figure represents corresponding values at each rotation of the impeller.

In this case one rotation needs one second of real time. It can be seen that, the

normalised variance computed based on mapping method deviates at the beginning

but afterwards approaches the CFD simulation values.

Fig. 5.11 shows the time evolution of normalised variance computed based on

mapping method and CFD simulation. Each marked point in this figure represents

corresponding values at each one third rotation of the impeller. Based on the geo-

metrical configuration of this particular kneader element we assume that the flow is

periodic after each one third rotation. It can be observed here that, the deviation

between the normalised variance computed based on mapping matrix and simulation,

is minimised. However, this deviation is caused by the different flow fields used for

mapping and tracking. The tracking data us created by introducing the particles at

t= 0 and then start the run. Where as the mapping data is created for the simulation

after the particles are uniformly distributed (after few impeller rotations) inside all

compartments. This means the flow field used for mapping is developed and the flow

field used for tracking is transient.

This is verified from the computations from Maxisch [52]. In Fig. 5.12 four data

sets Matrix1, Matrix2,... can be seen. Each data set represents the time evolution

of normalised variance for 1600 particles initially placed in one compartment. The

difference here is the creation of the matrix. All matrices are created for initially

randomly placed particles. The Matrix1 is created after the first rotation, Matrix2

is created after the second rotation and so on. As one can see, the ability of mixing

decreases when the flow field develops. In the simulation of the current thesis, the

mapping data is generated for a developed flow which according to the Fig. 5.12,

shows lower mixing effect. So the deviation between mapping and tracking (Fig. 5.10,

Fig. 5.11) can be seen. Better would be to compare tracking and mapping data for

either constant flow field or for transient flow field. In Fig. 5.13, the comparison of

tracking and mapping data was done for transient flow field. The agreement between

92 CHAPTER 5. MIXING IN A KNEADER ELEMENT

0.2

0.4

0.6

0.8

0 2 4 6 8

Rotation

Normalised varianc

Matrix1

Matrix2

Matrix3

Matrix4

Figure 5.12: Time evolution of normalised variance in a kneader element. Effect of

mapping matrices [52].

0.2

0.4

0.6

0.8

0 1 2 3 4 5Rotation

Normalised varianc

Tracking

Mapping

Figure 5.13: Time evolution of normalised variance in a kneader element. Comparison

of simulated data with data obtained from mapping matrix [52].

5.4. MAPPING METHOD 93

tracking and mapping data is very good. The evolution of normalised variance is

not decreasing as much as in the current thesis since the boundary conditions of the

simulations [52] are different in this particular case.

Thus, it can be concluded here that the mapping method is used to analyze the

mixture quality and optimize the mixers with less efforts since the matrix-vector multi-

plications only take a few seconds on a personal computer. Therefore, mapping method

can be used as an alternative engineering tool to conventional mixing measures.

94 CHAPTER 5. MIXING IN A KNEADER ELEMENT

Chapter 6

Conclusions and Outlook

Modeling and simulation of mixing in stirred vessels is very challenging task because

of complex design of impellers. It will be further complicated if the operation involves

highly viscous medium. The present contribution addressed the problems showed up

during the numerical simulation of highly viscous liquid mixing and recommended a

novel approach to overcome it.

Since this study deals with mixing in stirred vessels, various possible numerical

methods were discussed and validated by applying them to a simple rotating ref-

erence goemetry. The two well known, rotating reference frame and moving mesh

simulations were explained comprehensively and discussed the advantages and disad-

vantages. These two methods were employed on a model geometry and showed that

they both yield the same results. Further simulations were accomplished using the

rotating reference frame, because the efforts needed to set up rotating reference frame

are less compared with the moving mesh case set up and also moving mesh calculations

needs longer computational time. Additionally, a test case simulations were performed

in a flow between two concentric rotating cylinders (Couette-flow) to assess the capa-

bilities of Star-CD in handling highly viscous liquids. The torque computations for

highly viscous liquids in Coette-flow show remarkable agreement with the analytical

calculations.

Volume of Fluid (VOF) simulations in a partly filled anchor mixer reveals the flow

96 CHAPTER 6. CONCLUSIONS AND OUTLOOK

field distribution inside the mixing domain. The results show high velocity magnitudes

around the anchor arms as the flow is purely induced by the impeller. High pressure

in the front and low pressure behind the anchor arms attributes to the maximumu

forces acting on anchor arms against the flow. In addition, the systematic procedure

to calculate torque and power number in rotating systems is given with detailed math-

ematical expressions. Time evolution of torque in anchor mixer is shown for highly

viscous liquids Newtonian. The results show that torque is maximum at the onset

of motion due to the fluid at rest. However, it decreases abruptly and attains con-

stant value within a short period of time. Further, torque computations at different

viscosities also presented and discussed. It has to be noted that the power number

calculation also show the similar trend as torque, since power number is a function of

torque for a particular geometrical design.

Since the past two decades CFD has become a widely used tool for analysing, opti-

mising and supporting the design of mixing processes. There exists several ways for the

numerical investigation of mixing in stirred vessels. The well known among them are

numerical tracer experiments, Lagrangian particle tracking and entropy based mea-

sures. While numerical tracer experiments are of increasing importance as a means to

analyse mixing processes, there is a principle problem which is the so-called numerical

diffusion, i.e. the artificially generated smoothing of a tracer profile due to errors from

the discretization, and this numerical diffusion can be much stronger than the true

physical diffusion. The effect of numerical diffusion was accentuated with a proto type

lid driven cavity flow and highlighted the influence of solver settings and grid refine-

ment. A way to avoid this problem is to replace the continuous tracer concentration by

a number concentration obtained from Lagrangian (i.e. inertia free) particles that are

tracked during the simulation. This approach does not suffer from artificial diffusion,

since the position of tracer particles can be resolved with sub-grid-scale accuracy and

the velocity field at these particle positions can be obtained by interpolation from its

values at grid points.

In this investigation, the method of calculating intensity of segregation and hence

the intensity of mixing and mixing time are discussed based on the Lagrangian particle

tracking method. The total computational domain is divided into smaller compart-

ments (sub-volumes) and particles are initially placed in one compartment, say. During

the process of particle tracking the resulting number concentrations are recorded and

allow for computation of the evolution of its variance. A fundamental question in this

approach is how many compartments and particles are needed for a reliable assess-

ment of the mixing quality. Based on elementary statistics, it can be shown that a

reliable mixing time t1−τ

Mfor a given level  > 0requires 100/2particles (if standard

deviations instead of variances are employed) while a surprisingly small number of

about 20 compartments is sufficient.

This method has been evaluated using the numerical investigation of mixing in a

vessel stirred with an anchor impeller as well as a specific kneader element. It was

shown that even after a very long and expensive simulation, anchor mixer exhibits

poor mixing behaviour. It is known that anchor is very good for heat transfer between

vessel wall and fluid rather than its mixing capabilities. However, the simulations in a

specific kneader shows its excellent mixing characteristics within a short time. It is due

to the complex impeller design with pins on the kneading discs. It is believed that the

material inside the mixing chamber undergoes Baker transformation i.e., stretching

and folding.

Finally, we elaborated a mapping method to evaluate the quality of mixing. This

mapping method employs a transition matrix, which describes how many particles

are advected from one sub-volume to the other sub-volume in a particular period of

time. With the aid of this transition matrix one can compute variance evolutions and

mixing times using matrix-vector multiplications with significantly less computational

effort. Preliminary results of mapping method for a kneader element are discussed and

compared with the simulation results. The evolution of variance from both methods

(simulation and mapping) show good agreement between them. Therefore, the map-

ping method can be considered as an engineering tool to analyze the mixture quality

and optimize the mixers within few seconds on a personal computer. Nevertheless,

98 CHAPTER 6. CONCLUSIONS AND OUTLOOK

CFD simulation is performed to get the approximately equal particle distribution in

each compartments after one rotation.

In the present numerical investigations, it has been showed that Lagrangian par-

ticle tracking method is a novel approach to assess the mixing quality. The initial

position of blob of particles can be varied and checked for the influence on mixing

characteristics. Here, we restricted ourself to simple geometries. Generally, compu-

tation of mixing mechanisms in kneaders are difficult. Single shaft kneaders can be

calulated using appropriate moving or rotating reference frame as shown in this work.

But, the flow computation in double shaft kneaders cannot be modelled using the

sliding mesh technique because of the ovelapping of kneading discs as well as the close

clearance kneading pins. New techniques have been developed recently for the nu-

merical treatment of overlapping meshes. The above mentioned Lagrangian particle

tracking method can be extended to investigate mixing in double shaft kneaders and

also twin srew extruders.

Bibliography

[1] User Guide: Star-CD v 3.22, CD-Adapco (2004).

[2] Abid, M., Xuereb, C., and Bertrand, J. Hydrodynamics in vessels stirred

with anchors and gate agitators: Necessity of a 3-d modelling. Trans IChemE 70

(1992), 377–384.

[3] Anderson, F. B. U.S. Patent (1930).

[4] Anderson, J. D. Computational Fluid Dynamics: The basics with applications.

McGraw-Hill, Inc., 1995.

[5] Apruzzese, F., Pato, J., Balke, S. T., and Diosady, L. L. In-line mea-

surement of residence time distribution in a co-rotating twin-screw extruder. Food

Research International 36 (2003), 461–467.

[6] Aubin, J., and Xuereb, C. Design of multiple impeller stirred tanks for the

mixing of highly viscous fluids using cfd. Chem. Eng. Sci. 61 (2006), 2913–2920.

[7] Avalosse, T., and Rubin, Y. Analysis of mixing in corotating twin screw

extruders through numerical simulation. Int. Poly. Process. XV-2 (2000), 117–

123.

[8] Barailler, F., Heniche, M., and Tanguy, P. A. Cfd analysis of a rotor-

stator mixer with viscous fluids. Chem. Eng. Sci. 61 (2006), 2888–2894.

100 BIBLIOGRAPHY

[9] Beckner, J. L., and Smith, J. M. Anchor-agitated systems: Power input

with newtonian and pseudo-plastic fluids. Trans. Instn. Chem. Engrs. 44 (1966),

T224–T236.

[10] Bertrand, F., Tanguy, P. A., and Brito-De La Fuente, E. A new

perspective for the mixing of yield stress fluids with anchor impellers. Journal of

Chemical Engineering of Japan 29 (1996), 51–58.

[11] Bird, R. B., Stewart, W. E., and Lightfoot, E. N. Transport Phenomena.

Wiley, 1962.

[12] Boss, J. Evaluation of the homogeneity degree of a mixture. Bulk Solids Handling

6 6 (1986), 1207–1215.

[13] Bothe, D. Lecture Notes: Modelling of Chemical Reaction Systems. University

of Paderborn, WS 2004/05.

[14] Bothe, D., and Warnecke, H.-J. Berechnung und beurteilung strömungs-

basierter komplex-laminarer mischprozesse. Chem. Ing. Tech. 79 (2007), 1001–

1014.

[15] Carneiro, O. S., Covas, J. A., Ferreira, J. A., and Cerqueira, M. F.

On line monitoring of the residence time distribution along a kneading block of

a twin screw extruder. Polymer Testing 23 (2004), 925–937.

[16] Danckwerts, P. V. The definition and measurement of some characteristics

of mixtures. Appl. Sci. Res. A3 (1952), 279–296.

[17] de Graaf, R. A., Rohde, M., and Janssen, L. P. B. M. A novel model

predicting the residence time distribution during reactive extrusion. Chem. Eng.

Sci. 52 (1997), 4336–4345.

[18] Demirdzic, I., and Peric, M. Space conservation law in finite volume calcu-

lations of fluid flow. Int. J. Num. Meth. Fluids 8 (1988), 1037–1050.

BIBLIOGRAPHY 101

[19] Ebrahimi-Moshkabad, M., and Winterbottom, J. M. The behaviour of

an intermeshing twin screw extruder with catalyst immobilised screws as a three

phase reactor. Catalysis Today 48 (1999), 347–355.

[20] Elemans, P. H. M., and Meĳer, H. E. On the modeling of continuous

mixers. part ii: The cokneader. Poly. Eng. Sci. 30 (1990), 893–904.

[21] Ferziger, J. H., and Peric, M. Computational Methods for Fluid Dynamics.

Springer, 1999.

[22] Fogiel, M. The transport phenomena problem solver: Momentum, Mass and

Energy. Research and Education Association, 1986.

[23] Galaktionov, A. S., Anderson, A. D., and Peters, G. W. M. Mixing

simulations: tracking strongly deforming fluid volumes in 3d flows. In M Bubak

JD, Wasniewski J (eds) Lecture notes in computer science, volume 1332 of Recent

advances in parallel virtual machine and message passing interface. Springer,

Heidelberg. 463-469 (1997).

[24] Galaktionov, O. S., Anderson, P. D., Kruĳt, P. G. M., Peters, G.

W. M., and Meĳer, H. E. M. A mapping approach for three-dimensional

distributive mixing analysis. Computers & Fluids 30 (2001), 271–289.

[25] Gaspar-Cunha, A., Poulesquen, A., Vergnes, B., and Covas, J. A.

Optimization of processing conditions for polymer twin-screw extrusion. Int.

Poly. Process. XVII-3 (2002), 201–213.

[26] Giguere, R., Bertrand, F., and Tanguy, P. A. A three-dimensional mesh

refinement strategy for the simulation of fluid flow with fictitious domain method.

Comp. & Chem. Eng. 30 (2006), 453–466.

[27] Glowinski, R., Pan, T., and Periaux, J. A fictitious domain method for

dirichlet problem and applications. Comp. Meth. Appl. Mech. Eng. 111 (1994),

283–303.

102 BIBLIOGRAPHY

[28] Glowinski, R., Pan, T.-W., Helsa, T., and Joseph, D. A distributed

lagrange multiplier/fictitious domain method for particulate flows. Int. J. Multi-

phase Flows 25 25 (1999), 755–794.

[29] Grebe, T. Simulation und Modellierung des Mischverhaltens von Taylor-Couette

Reaktoren. PhD thesis, University of Paderborn, 2004.

[30] Hadzic, H. Development and Application of a Finite Volume Method for

the Computation of Flows Around Moving Bodies on Unstructured, Overlapping

Grids. PhD thesis, Technical University of Hamburg-Harburg, 2005.

[31] Harnby, N., Edwards, M. F., and Nienow, A. W. Mixing in the Process

Industries. Butterworth-Heinemann, 1992.

[32] Iranshahi, A., Heniche, M., Bertrand, F., and Tanguy, P. A. Numer-

ical investigation of the mixing efficiency of the ekato paravisc impeller. Chem.

Eng. Sci. 61 (2006), 2609–2617.

[33] Kaminoyama, M., Arai, K., and Kamiwano, M. Numerical analysis of

power consumption and mixing time for a pseudoplastic liquid in geometrically

similar stirred vessels with several kinds of plate-type impellers. Journal of Chem-

ical Engineering of Japan 27 (1994), 17–24.

[34] Kaminoyama, M., Saito, F., and Kamiwano, M. Flow analogy of pseudo-

plastic liquid in geometrically similar stirred vessels based on numerical analysis.

Journal of Chemical Engineering of Japan 23 (1990), 214–221.

[35] Kang, T. G., and Kwon, T. H. Colored particle tracking method for mixing

analysis of chaotic mixers. Journal of Micomech. Microeng. 14 (2004), 891–899.

[36] Kresta, S. M., Krebs, R., and Martin, T. Review: The future of mixing

research. Chem. Eng. Tech. 27 (2004), 208–214.

BIBLIOGRAPHY 103

[37] Kruĳt, P. G. M. Analysis and optimization of laminar mixing: design, develop-

ment and application of the mapping method. PhD thesis, Eindhoven Universtiy

of Technology, 2000.

[38] Kruĳt, P. G. M., Galaktionov, O. S., Anderson, P. D., Peters, G.

W. M., and Meĳer, H. E. M. Analyising fluid mixing in periodic flows by

distribution matrices: the mapping method. AIChE Journal 47 (2001), 1005–

1015.

[39] Kuriyama, M., Inomata, H., Arai, K., and Saito, S. Numerical solution

for the flow of highly viscous fluid in agitated vessel with anchor impeller. AIChE

Journal 28 (1982), 385–391.

[40] Lacely, P. M. C. Development in the theory of particle mixing. Journal of

Applied Chemistry 4 (1954), 257–268.

[41] Leonard, B. P. A stable and accurate convective modelling procedure based

on quadratic upstream interpolation. Comp. Meth. Appl. Mech. Eng. 19 (1979),

59–98.

[42] Lertwimolnun, W., and Vergnes, B. Effect of processing conditions on the

formation of polypropylene/organoclay nanocomposites in a twin screw extruder.

Poly. Eng. Sci. 46 (2006), 314–323.

[43] List, H. Swiss Patent (1945).

[44] Liu, W., and Haller, G. Strange eigenmodes and decay of variance in the

mixing of diffusive tracers. Physica D 188 (2004), 1–39.

[45] Lyu, M. Y., and White, J. L. Models of flow and experimental studies on a

modular list/buss kneader. Int. Poly. Process. 10 (1995), 305–313.

[46] Lyu, M. Y., and White, J. L. Modeling of a viscous non-newtonian polymer

melt in a list/buss kneader and comparison to experiment. Int. Poly. Process. 11

(1996), 208–221.

104 BIBLIOGRAPHY

[47] Lyu, M. Y., and White, J. L. Simulation of linear viscoelastic flow behavior

in the buss kneader. Poly. Eng. Sci. 37 (1997), 623–635.

[48] Lyu, M. Y., and White, J. L. Simulation of non-isothermal flow in a modular

buss kneader and comparison with experiment. Int. Poly. Process. 12 (1997),

104–109.

[49] Lyu, M. Y., and White, J. L. Residence time distribution and basic studies

of flow and melting in a modular buss kneader. Poly. Eng. Sci. 38 (1998), 1366–

1377.

[50] Mackley, M. R., and Saraiva, R. M. C. N. The quantitative description

of fluid mixing using lagrangian- and concentration-based numerical approaches.

Chem. Eng. Sci. 54 (1999), 159–170.

[51] Maridass, B., and Gupta, B. R. Performance optimization of a counter rotat-

ing twin screw extruder for recycling natural rubber vulcanizates using response

surface methodology. Polymer Testing 23 (2004), 377–385.

[52] Maxisch, T. Private communication. Tech. rep., Institute of Technical Chem-

istry and Chemical Process Engineerng, University of Paderborn, 2008.

[53] Mehranpour, M., Nazokdast, H., and Dabir, B. Computational study of

the velocity field in the conveying element of a ko-kneader with cfd method. Int.

Poly. Process. 17 (2002), 108–114.

[54] Meĳer, H. E., and Elemans, P. H. M. The modeling of continuous mixers.

part i: The corotating twin-screw extruder. Poly. Eng. Sci. 28 (1988), 275–290.

[55] Mezaki, R., Mochizuki, M., and Ogawa, K. Engineering data on mixing.

Elsevier, 2000.

[56] Motzigemba, M., Bothe, D., Broecker, H.-C., Prüss, J., and War-

necke, H.-J. A contribution to simulation of mixing in screw extruders em-

BIBLIOGRAPHY 105

ploying commercial cfd-software. In 10th European conference on mixing, Delft

(2000).

[57] Muralidhar, K., and Biswas, G. Advanced Engineering Fluid Mechanics.

Narosa Publishing House, 1999.

[58] Ottino, J. M. The kinematics of mixing: streching, Chaos, and Transport.

Cambridge University Press, Cambridge, 1989.

[59] Pedrosa, S. M. C. P., and Nunhez, J. R. The behavior of stirred vessels

with anchor type impellers. Comp. & Chem. Eng. 24 (2000), 1745–1751.

[60] Peters, D. C., and Smith, J. M. Fluid flow in the region of anchor agitated

blades. Trans. Instn. Chem. Engrs. 45 (1967), T360–T366.

[61] Phelps, J., and Tucker III, C. L. Lagrangian particle calculations of dis-

tributive mixing: Limitations and applications. Chem. Eng. Sci. 61 (2006), 6826–

6836.

[62] Prat, L., Guiraud, P., Rigal, L., and Courdon, C. Two phase residence

time distribution in a modified twin screw extruder. Chem. Eng. Processing 38

(1999), 73–83.

[63] Puaux, J. P., Bozga, G., and Ainser, A. Residence time distribution in a

corotating twin-screw extruder. Chem. Eng. Sci. 55 (2000), 1641–1651.

[64] Schütz, S., Bierdal, M., and Piesche, M. Charakterisierung des mis-

chverhaltens von gegenstrom-injektions-mischern. Chem. Ing. Tech. 77 (2005),

398–405.

[65] Seck, O., and Maxisch, T. Verfahrenstechnische charakterisierung eines

kneters. Master’s thesis, University of Paderborn, 2006.

[66] Seidel, W. Note on a metrically transitive system. In Proceedings of National

Academy of Science USA 19 (1933).

106 BIBLIOGRAPHY

[67] Shon, K., Chang, D., and White, J. L. A comparative study of residence

time distribution in a kneader, continuous mixer, and modular intermeshing co-

rotating and counter rotating twin screw extruders. Int. Poly. Process. 14 (1999),

44–50.

[68] Singh, M. K., Kang, T. G., Meĳer, H. E. H., and Anderson, P. D.

The mapping method as a toolbox to analyze, design, and optimize micromixers.

Microfluidics and Nanofluidics DOI 10.1007/s10404-007-0251-7.

[69] Sombatsompop, N., and Panapoy, M. Comments on temperature profiles

of pp melt in the barrel of a twin screw extruder. Polymer Testing 20 (2000),

217–221.

[70] Sommerfeld, M., and Decker, S. State of the art and future trends in cfd

simulation of stirred vessel hydrodynamics. Chem. Eng. Tech. 27 (2004), 215–224.

[71] Spencer, R. S., and Wiley, R. M. The mixing of very viscous liquids. Journal

of Colloid Science 6 (1951), 133–145.

[72] Spurk, J. H., and Aksel, N. Strömungslehre. Springer, 2005.

[73] Stuber, N. P., and Tirrell, M. Continuous polymerization studies in a twin

screw extruder. Poly. Proc. Eng. 3 (1985), 71–83.

[74] Sumi, Y., and Kaminawo, M. Development and mixing characteristics of

a multistage impeller for agitating highly viscous fluids. Journal of Chemical

Engineering of Japan 34 (2001), 485–492.

[75] Thiffeault, J.-L. Scalar decay in chaotic mixing. transport in geophysical

flows: Ten years after. In Proceedings of the Grand Combin Summer School

(2004).

[76] Trolestra, E. J., Van Dierendronck, L. L., and Janssen, L. P. B. M.

Modeling of a buss-kneader as a polymerization reactor for acrylate. parti: Model

validation. Poly. Eng. Sci. 42 (2002), 230–239.

BIBLIOGRAPHY 107

[77] Trolestra, E. J., Van Dierendronck, L. L., and Janssen, L. P. B. M.

Modeling of a buss-kneader as a polymerization reactor for acrylate. partii:

Methyl methacrylate based resins. Poly. Eng. Sci. 42 (2002), 240–247.

[78] Unlu, E., and Faller, J. F. in twin screw food extrusion. Journal of Food

Engineering 53 (2002), 115–131.

[79] van Zuilichem, D. J., Kuiper, E., Stolp, W., and Jager, T. Mixing ef-

fects of constituting elements of mixing screws in single and twin screw extruders.

Powder Technology 106 (1999), 147–159.

[80] Vanyo, J. P. Rotating Fluids in Engineering and Sciences. Butteworth-

Heinemann, 1993.

[81] Versteeg, H. K., and Malalasekara, W. An Introduction to Computa-

tional Fluid Dynamics. Longman Scientific & Technical, 1995.

[82] Wang, W., Manas-Zloczower, I., and Kaufman, M. Entropic character-

isation of distributive mixing in polymer processing equipment. AIChE Journal

49 (2003), 1637–1644.

[83] White, J. L., and Potente, H. Screw Extrusion. Carl Hanser Verlag, 2003.

[84] Wurster, C. German Patent (1901).