Three-Dimensional Flow Analysis in The Calendering Process

Download as pdf or txt
Download as pdf or txt
You are on page 1of 92

University of Thessaly

School of Engineering
Department of Mechanical Engineering

Three–Dimensional Flow Analysis in the


Calendering Process

A Thesis
Submitted in partial fulfillment of the requirements for the degree of
Master of Science (M.Sc.)
Of the Department of Mechanical Engineering of University of Thessaly

by

Nickolaos D. Polychronopoulos
B.Sc. in Materials Science

Advisory Committee:
Assoc. Prof. T.D. Papathanasiou, supervisor
Prof. N. Pelekasis
Prof. P. Tsiakaras

Volos 2012

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
© 2012 Πολυχρονόπουλος Νικόλαος

Η έγκριση της μεταπτυχιακής εργασίας από το Τμήμα Μηχανολόγων Μηχανικών της


Πολυτεχνικής Σχολής του Πανεπιστημίου Θεσσαλίας δεν υποδηλώνει αποδοχή των
απόψεων του συγγραφέα (Ν. 5343/32 αρ. 202 παρ. 2).

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Τὰ πάντα ῥεῖ, τὰ πάντα χωρεῖ καὶ οὐδὲν µένει.

Ηράκλειτος, 554 π.Χ. – 484 π.Χ.

ii

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Acknowledgements

I am most grateful to my supervisor Assoc. Prof. Thanasis D. Papathanasiou for


his guidance, criticism and continuous support during my Master’s Degree. I would like
to thank Prof. Nikolaos Pelekasis, for his valuable knowledge in Fluid Mechanics and
computational algorithms for fluid flows. I would also like to thank Prof. Panagiotis
Tsiakaras and Assoc. Prof. Nikolaos Andritsos for their crucial contribution in
understanding the true meaning of physical phenomena. I would like to thank Prof.
Emeritus of McMaster University, John Vlachopoulos, for his guidance and continuous
effort to conceive in the best possible way on how nature works and getting me one step
further, beyond the mathematical equations. Very special thanks go to Dr. Ioannis Sarris,
Dr. Sotirios Kakarantzas and Dr. Christos Dritselis for their help in understanding the
meaning of Computational Fluid Dynamics. I would like also to thank my family for their
true support and my friends, Petros Christodoulou and Dimitris Tsokolis for their
memorable company and friendship. Last but not least, I would like to thank, the PhD
candidate and my friend, Lefteris Benos, for the inexhaustible conversations to
understand the strange paths of nature.

iii

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Objectives

The importance to study the calendering process lies to the fact of the high
volume capacities and that the formed sheet can have very high quality. Thus, operating
variables used in engineering calculations need to be evaluated beforehand with the most
important being the maximum pressure that develops in the process. Pressure distribution
is developed in the machine and cross–machine direction, primarily due to the rotational
movement of the calenders. A lot of work has been carried out in the past on an effort to
study the flow numerically using one–dimensional or two–dimensional models as well as
experimentally. However, there are very few studies concerning the three–dimensional
effects occurring during the physical process. These effects, which in the past were
observed experimentally, are the pressure distribution in the cross–machine direction
which is described by a quasi–parabolic pressure profile, the spiraling flow pattern of the
material in the melt bank and the spreading of the sheet as it passes through the calenders.
The scope and the contribution of the present work is to employ the Generalized
Newtonian Fluid model and solve the fully 3D Navier–Stokes equations so that to predict
the three–dimensional effects as stated above. This was carried out by making several
assumptions and approximations concerning the rheological behavior of the material as
well as the computational flow domain. All the numerical experiments were carried out
using the Open Source Field Operation and Manipulation (OpenFOAM) software.

iv

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
THESIS OUTLINE

In Chapter 1 a brief overview concerning the calendering process is presented,


explaining the way polymeric melts form thin sheets by the process. Chapter 2 presents
the mathematical formulation used to describe the process by means of the Lubrication
Approximation Theory (LAT) and how the equations under the LAT can be derived from
the well-known Navier–Stokes equations making assumptions concerning the physical
behavior of polymer melts. Under certain flow conditions and geometries the Navier–
Stokes can be simplified. The formulation is done for Newtonian and non–Newtonian
fluids. Chapter 3 presents a literature survey concerning previous studies on the process,
with the pressure distribution in the machine direction and the pattern of the flow being
the most important aspects studied. The survey covers one-dimensional and two-
dimensional analyses of the calendering process. In Chapter 4, the Finite Volume Method
that OpenFOAM employs, is briefly discussed. It is presented the type of equations the
software solves and how these equations were solved for the study of the calendering. In
Chapter 5, the analytical solution obtained from the Lubrication Approximation Theory is
compared against the numerical solution of the Navier–Stokes equations, using the open
source software OpenFOAM. We consider both, Newtonian and Generalized–Newtonian
Fluids utilizing the power–law and Carreau viscosity models. The dimensions of the
calendering apparatus as well as the rheological parameters of the fluids were obtained
from the literature and represent actual industrial equipment and real materials. The
polymeric material was Poly–VinylChloride (PVC). Finally, in Chapter 6, we focus on
the three-dimensional flow features of the process. There is only a limited amount of
published work concerning 3D calendering and this is presented and discussed here, in
light of our own 3D computational results for Newtonian and non–Newtonian fluids.

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Table of Contents

Chapter 1 ....................................................................................................... 1
Introduction ..................................................................................................................... 1
Chapter 2 ....................................................................................................... 4
Mathematical Formulation of Calendering ..................................................................... 4
2.1 General Considerations ......................................................................................... 4
2.2 Lubrication Approximation with the Newtonian Model .................................... 11
2.3 Lubrication Approximation with the power–law Model .................................... 16
Chapter 3 ..................................................................................................... 19
Literature Survey of Calendering.................................................................................. 19
Chapter 4 ..................................................................................................... 25
Finite Volume Method .................................................................................................. 25
Introduction ............................................................................................................... 25
4.1 Discretisation of Continuous Operations ............................................................ 25
4.2 Governing and Constitutive Equations of the Solver ......................................... 27
4.3 Discretised Equations.......................................................................................... 28
4.3.1 System of Equations .................................................................................... 28
4.3.2 Linearization of the Navier – Stokes Equations .......................................... 30
4.4 Numerical Setup.................................................................................................. 32
Chapter 5 ..................................................................................................... 33
Calendering Validation with OpenFOAM .................................................................... 33
Introduction ................................................................................................................... 33
5.1 Newtonian Fluid.................................................................................................. 33
5.1.1 Flow Domain Geometry .............................................................................. 33
5.1.2 Boundary Conditions ................................................................................... 36
5.1.3 Results and Discussion ................................................................................ 38
5.2 Generalized Newtonian Fluid Model .................................................................. 47
5.2.1 Flow Domain Geometry and Boundary Conditions .................................... 47
5.2.2 Results and Discussion ................................................................................ 48
Chapter 6 ..................................................................................................... 53
3D Effects in Calendering ............................................................................................. 53
Introduction ............................................................................................................... 53

vi

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
6.1 Numerical Results for a Newtonian Fluid .......................................................... 58
6.1.1 Geometry of Computational Flow Domain and Boundary Conditions ....... 58
6.1.2 Results and Discussion ................................................................................ 62
6.2Numerical Results with Carreau Viscosity Model............................................... 71
6.2.1 Geometry, Spreading of the Sheet and Boundary Conditions ..................... 71
6.2.2 Results and Discussion ................................................................................ 72
Closure and Future Work .......................................................................... 79

vii

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Chapter 1

Introduction

The term “calender” is derived from the Greek word Κύλινδρος (cylinder) and
according to Merriam–Webster’s International Dictionary, it means “to press (as cloth,
rubber or paper) between rollers or plates in order to smooth and glaze or to thin into
sheets”. Calendering is the process in which molten material is dragged through the
narrow region between two rotating (counter rotating or co–rotating) heated rolls in such
a way as to produce a film or sheet. Figure 1 shows a schematic representation of the
calendering process. It is observed that due to the reduced area, a melt bank is created
before the region where the distance between the rollers is the smallest (the nip region).
In this melt bank, a very interesting flow pattern develops with multiple recirculation
regions (this is shown later on the thesis). The number of rolls of a calendering line is
determined primarily by the nature of the material processed as well as the desired
product. For example, rubber can be processed and calendered on a two–roll calendar,
with four roll calenders generally used for double coating of substrates (Figure 1.2.a).
However, the surface requirements for good to excellent quality of calendered
thermoplastic polymer require four–roll calenders as shown at Figures 1.2.b and 1.2.c
which are the standard calender configurations typically called “inverted L” or “Z”
configurations. Calender sizes range up to 900 mm in diameter and 2500 mm wide, with
polymer throughputs up to 4000 kg/h. The surface temperature of the rollers is carefully
controlled by using drilled rolls (i.e. axially drilled holes all around the periphery) in
which a temperature – controlling liquid is circulated [1,2,3].

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 1.1. Sketch showing the calendering process between two counter–rotating rolls [3].

In an analysis of the calendering process one seeks the relationships among


performance variables, such as the sheet thickness as well as design and operational
variables, such as roll diameter and the speed of the rollers. Also, the effect of the
different rheological parameters on these relationships is sought. Calenders may also be
used for the production of a certain “finish” on a surface, either a degree of “gloss” or
smoothness which might be required for a desired use, or even produce a sheet with
intentional degree of roughness or pattern. Finish is a result of interaction between the
fluid and the surface of the rollers in the region where the calendered sheet separates
from the rollers [1].

Figure 1.2. (a) A four–roll, inclined “Z” calender for double casting of tire cord. (b) A four–roll, inverted
“L” calender. (c) A four – roll “Z” calender [2].

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
The calendering process, is commonly used for shaping high melt viscosity
thermoplastic sheets and is particularly suitable for polymers susceptible to thermal
degradation or containing substantial amounts of solid additives. This is due to the fact
that the calender can convey large rates of melt with a small mechanical energy input (in
comparison to a single–screw extruder for example).

A very important aspect in the calendering process is the thickness uniformity of


the produced sheet in both the machine and cross – machine direction. Variations in gap
size due to the dimensions of the roll, thermal effects, and roll distortions due to high
pressures developing in the gap, will result in product nonuniformity in the cross–
machine direction. Eccentricity of the roll with respect to the roll shaft, as well as roll
vibration and feed uniformity, must be tightly controlled to avoid non - uniformity in the
machine direction. A uniform empty gap size will be distorted in operation because of
hydrodynamic forces, developed in the nip region, which deflect the rolls. From such a
condition, the resulting product will be thick in the middle and thin at the edges.

Calendering lines are very expensive in terms of capital investment in machinery.


Film and sheet extrusion are competitive processes because the capital investment for the
extruder itself is only a fraction of the cost of a calender. However, the high quality and
volume capabilities of calendering lines make them advantageous for many types of
products, especially for temperature sensitive materials. Polyvinylchloride (PVC) is the
major polymer that is calendered.

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Chapter 2

Mathematical Formulation of Calendering


2.1 General Considerations

Polymer melts flow at very low Reynolds numbers (e.g. Re = 10-2–10-4) and the
creeping flow approximation is applicable, in which, the intertial forces are very small
compared to the viscous forces and the flow is steady in time. Thus, the conservation
equations can be simplified into the following

  U  0 (2.1)

0  p    (2.2)

C p U  T    kT     U   : U (2.3)

where ρ is the density, U = U(ux,uy,uz) the velocity vector, p the pressure,  the stress
tensor, Cp specific heat, T temperature and k thermal conductivity. The last commonly
used set of approximations in liquid flow (isothermal, as well as nonisothermal) and in
conductive heat transfer, is to treat k, Cp, and ρ as constant quantities, independent of T
and p. In polymer processing, where both heat transfer and flow take place, typical
pressure variations may reach up to 40 MPa for most extrusion processes, such as flat die
extrusion and this suggests density and viscosity changes. While the effect of pressure in
melt density is small (perhaps 5% under the most severe extrusion conditions) the effect
on viscosity is larger. Cogswell [4] suggests as a very rough approximation that 10 MPa
of pressure have as much effect as a reduction of 5oC of temperature. While, the viscosity

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
dependence on temperature is nearly always taken into consideration the author is not
aware of any serious attempts to consider pressure dependence of viscosity in the
extrusion and the calendering process. For those reasons, in the presentation below the
incompressibility assumption will be used and the viscosity will be a function of shear
rate but not of pressure. In the Generalized Newtonian Fluid (GNF) model, which will be
presented later on the present study, the viscosity is always a function of temperature,
usually through an exponential equation. However, in the present study any thermal
effects of the process were neglected, and it was decided that the viscosity will be only a
function of the shear rate. Thus, Eq. 2.1, 2.2 and 2.3 can be written in the following form

U  0 (2.4)

0  p    (2.5)

C p U  T  k2T   : U (2.6)

Moreover, assuming that the fluid does not spread as it enters the gap between the rolls,
we may write Eq. 2.4, 2.5 and 2.6 in two dimensions, where x is the direction of the flow
and y perpendicular to roll axis (thus neglecting phenomena in the direction of the
cylinder axis), as

u x u y
 0 (2.7)
x y

p  xx  xy
0   (2.8)
x x y

p  yx  yy
0   (2.9)
y x y

 T T    2T  2T  u u u u
C p  ux  uy   k  2  2    xx x   xy x   yx y   yy y
 x y   x y  x y x y (2.10)

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Furthermore, it can be assumed that the most important effect in the calendering process,
namely the build-up of pressure and the orientation of the polymeric molecules, occur in
the region of the minimal roll separation, which is called nip region.

centerline

Figure 2.1. Schematic representation for definition of the calendering flow analysis with the Lubrication
Approximation Theory.

In that region and extending to either side (in the ±x direction) by a distance of the order
of xf as shown in Figure 2.1, the roll surfaces are nearly parallel if, as is usually the case
Ho << R. For that reason, it is reasonable to assume that the flow in the gap will be nearly
parallel so that the key assumptions of the Lubrication Approximation Theory (LAT) may
be applied

 
 , u y  u x , u x  u x x, y  and p  p(x)
x y

thus we have

u x u y
 0 (2.11)
x y

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
dp  xy
0  (2.12)
dx y

T  2T u
C pux  k 2   xy x
x y y (2.13)

The process of transforming problems from Eq. 2.8 and 2.9 into Eq. 2.12 is often referred
to as Lubrication Approximation. The reference to lubrication comes from the fact that
lubrication problems themselves typically involve a geometry such that Eq. 2.12 is valid.
The continuity equation (Eq. 2.11) may be replaced by the integral form

h( x)

Q  2  u x dy (2.14)
0

where Q is the volumetric flow rate.

For the above equations to be solved a constitutive equation that relates stress to
the rate of strain is also needed. Constitutive equations are relations between stresses and
strains (deformations) and are solved together along with the continuity and momentum
equations. In its simplest form (simple shear flow), the Newtonian equation is a linear
relation between the shear stress and the shear rate

du x
     (2.15)
dy

where  is the shear rate. This expression for  is valid for simple shear flows between
two flat plates for instance, and it is directly applicable to unidirectional flows. In
polymer processing, however, most interesting flow problems require two– or three–
dimensional analyses (like the present study), of creeping (low Reynolds, Re<<1) flows
and the shear rate is expressed in a tensor form that will be subsequently presented.

The stress tensor  , in Eq. 2.5 is a (second order) tensor as shown below

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
 xx  xy  xz 
 
   yx  yy  yz  (2.16)
 
 zx  zy  zz 

To generalize the Newtonian equation in three dimensions we must propose a


linear relation between stress components and rate of strain components. The rate of
strain tensor is the following

1  u u 
D
1
2
 
u  uT   i  j 
2  x j xi 
(2.17)

Where the components of the velocity gradient tensor, u are the following

ui
ui , j  (2.18)
x j

and

 Dxx Dxy Dxz 


 
D   Dyx Dyy Dyz  (2.19)
D Dzy Dzz 
 zx

where

1  u u  1  u u 
Dxx   x  x  , Dxy   x  y  (2.20)
2  x x  2  y x 

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
and similarly the other components can be written out explicitly in terms of the
components in the x, y and z directions. The Newtonian constitutive equation can be
generalized in the form

    2 D  (2.21)
 

and this means that

ux  u u y 
 xx   2Dxx   2 ,  xy   2 Dxy     x   etc. (2.22)
x  y x 

where in the above equations η is the fluid viscosity, which, under the assumption of a
Newtonian fluid is constant, independent on both, shear rate and temperature. To account
for the shear–thinning behavior of polymer melts, a number of models have been
proposed in the literature. These include the power–law, Carreau–Yasuda (the model will
be presented further on the present study) and Cross models, are generalized by replacing

 by the function (IID) which is the second invariant of the strain rate tensor 2 D . Thus
we have the power–law equation which can be written as

  mT   
n 1
(2.23)

where η is the viscosity, m is the consistency index (the larger the m the more viscous the
polymer melt), n is the power–law index, expressing the strength of the influence of shear
rate on viscosity (for n=1 we have Newtonian behavior) and  is the magnitude of the

rate of strain tensor and is given by the following equation

1
  ΙΙ D (2.24)
2

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
A severe limitation of the power-law model and the source of large errors when
used in 2D or 3D flows is that at very low shear rates, it predicts abnormally large
viscosities and thus flow fields which can be severely unrealistic.

Other viscosity models which are widely used and do not share the above
limitation, are the Carreau – Yasuda and Cross models. The Carreau – Yasuda model is
the following

     1    a 


n 1
(2.25)
 

For α=2 the model is simply called Carreau. The Cross model is the following


    (2.26)
1   
1n

where η0 is the zero–shear rate viscosity. The parameter λ is a constant with units of time.
In the Carreau–Yasuda model, λ determines the shear rate at which a transition occurs
from the zero–shear rate plateau to the shear thinning portion of the viscosity versus shear
rate curve. In the present study, except the Newtonian model, power–law and Carreau
viscosity models will be used.

In the power–law model the consistency index m (i.e. the value of the viscosity at
shear rate   1s 1 ) is a function of temperature in the same way as zero–shear viscosity.
The power–law n parameter is usually not a function of temperature. The viscosity
dependence of temperature is usually expressed either by an Arrhenius expression,
mostly used in polymer physics and rheology

E  1 1 
   ref .     (2.27)
 R  T Tref . 

10

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
or simple exponential, mostly used in equipment design

  ref exp[bT  Tref . ]


(2.28)

where ηref. is a viscosity measured at a reference temperature, E is the activation energy


and R is the gas constant. The temperature sensitivity coefficient b is usually between
0.01 and 0.1 oC-1 for most commercial polymers. For HDPE (linear polymer) the value of
b is roughly 0.01, while for LDPE (branched) it may reach 0.03.

2.2 Lubrication Approximation with the Newtonian Model

For a Newtonian fluid the stress is related to the rate of strain in the following
way

u x
 xy   (2.29)
y

where μ is the Newtonian viscosity. The momentum equation with the Lubrication
Approximation analysis for the rectangular coordinate system shown at Figure 2.1 is

p  2u
  2x (2.30)
x y

The above equation can be solved with the no–slip and symmetry boundary conditions

u x
u x  U on y  hx  and  0 on y  0 (2.31)
y

where U is the velocity of the calenders, and they give

11

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
y 2  h 2 x  dp
ux  U  (2.32)
2 dx

Examining h(x), which is the y distance from the centerline, it is easy to verify the
change of that distance with respect to the position x along the centerline by the following
expression


hx   H o  R  R 2  x 2 
12
(2.33)

However, in order to simplify the analysis the above equation can be expressed with a
remarkable degree of accuracy by a second order polynomial since it is more likely to
confine the analysis to values that x<<R. Thus Eq. 2.33 can be approximated as

 x2 
hx   H  1   (2.34)
 2 H o R 

The above equation actually treats the surfaces of the calenders as parabolic. It is
convenient now to define a set of dimensionless variables as below

x y
x' y' 
2 RH o Ho

ux pH o
ux '  p' 
U U

Then, we can find an expression for the pressure gradient from the integrated form of the
continuity equation (Eq. 2.14) like below

 h 2 x  dp 
h( x)

Q  2  u x dy  2hx U   (2.35)


0  3 dx 

12

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Then the equation above can be solved for dp/dx and find in dimensionless form

dp' 18R x'2 2


 (2.36)
dx' 
H o 1  x' 2 3 

where λ is a dimensionless parameter

H ( x)
2  1 (2.37)
Ho

From a physical point of view, the effects occurring at the detachment point (i.e. at y=Hf
and x=xf), will exert a strong influence on the overall behavior of the process. Eq. 2.36
gives zero pressure gradient at x'   . The distance x'   represents the point where the
sheet leaves the rolls (the point of detachment) and x'   is the point of maximum
pressure. The point of detachment is obtained if we ignore the forces acting at the
separation boundary and then we expect the stress just inside the fluid to be the same as
the ambient pressure. It is also assumed that right at separation, the profile is flat and the
viscous stresses vanish at that point. The flat profile corresponds to u x y  0 and
2
 2u x y  0 . For those reasons the boundary condition for the pressure at the

detachment point is the following

p' 0 at x'   (2.38)

Eq. 2.36 may be integrated with respect to x' and if the boundary condition of Eq. 2.38 is
imposed, then the solution for p' x' is found to be the following

p' 
9R

 
 x'2 1  32  1  52
x ' 
1  32

tan 1
x ' tan 
1
 
1  32 
 (2.39)
32 H o  1  x' 22
1  2 

13

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 2.2, shows qualitatively the shape of the pressure distribution along the centerline.
The maximum pressure occurs just upstream, at distance x'   from the nip region and

Figure 2.2. Schematic representation of the pressure distribution between the calenders as sketched by
Middleman [1].

has the value of

3C   R
pmax '  (2.40)
2 2H o

where

1  32
C   
1  2
 
  1  32 tan 1  (2.41)

and becomes zero at x'   where the polymer detaches from the calenders. At that point
the pressure gradient vanishes. There is also another value of x' at which the pressure
gradient vanishes and is seen to be at x'   . The pressure there is obtained if we set
x'   in Eq. 2.39 from which it is obtained

1  32 2  1 
p' x   
9R

 
  1 3   tan   (2.42)
 
2
32 H o  1 2

14

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
The most reasonable assumption to be made here is

p' x    0 (2.43)

from which is followed that the dimensionless parameter λ has a finite value, namely

o  0.475 (2.44)

and substituting that value in Eq. 2.43 then it can be seen that

Hf
 1    1.226
2
(2.45)
Ho

and that the sheet thickness depends only on Ho. At that point it is noted that λο is the
dimensionless leave–off distance when Eq. 2.43 is imposed. Thus, using the Lubrication
Approximation Theory the final sheet thickness can be evaluated by measuring the
distance Ho of the calenders at the nip region. The velocity and the pressure profiles can
be derived with the assumption of infinite reservoir far upstream behind the nip region.

Figure 2.3. Recirculation flow pattern in the nip region. In the schematic only the half–region is shown
from the top roller to the centerline [1].

Also the velocity profile can be easily obtained from Eq. 2.32 with the help of Eq. 2.36
and can be expressed as

15

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
3 1  y'2 2  x'2 
ux  1
2 1  x'2  (2.46)

The velocity distribution across the gap through the calenders can be seen at Figure 2.3.
Because of the pressure build-up at the nip region, there is a backflow component for
x'   which is superimposed onto the drag flow component. Thus, near the entrance of
the region between the calenders a recirculatory flow pattern develops. More on this topic
will be discussed in a subsequent chapter.

2.3 Lubrication Approximation with the Power–Law Model

For the power–law viscosity model, the derivation of the pressure solution
according to the LAT proceeds in a manner similar to that outlined earlier for the
Newtonian fluid. For a power–law fluid, the shear stress τxy can be expressed in the
following simple form

n 1
u u x
 xy m x (2.47)
y y

where the absolute value sign on the velocity gradient term avoids the problem of taking
a fractional root of a negative number.

From the Lubrication Approximation Theory solution for the Newtonian model we
anticipate two flow regions, a region where the pressure gradient is positive    x'  
and a region where the pressure gradient is negative    x'   . Thus, integrating Eq.
2.12 by means of Eq. 2.47 separately for each region we have

1
1  1 dp  n  1 n n 
n
ux  U      y  hx 1 n  (2.48)
q  m dx   

16

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Introducing either of the above equations into the integrated equation of the conservation
of mass (Eq. 2.14) and then solving for the pressure gradient, we have

dp'  2n  1 
n

2 R   x'   x'
2 2 2
 2 n 1

   (2.49)
dx'  n  Ho 
1  x' 2
2 n 1

where the dimensionless pressure is the following

n
pH 
p'   o  (2.50)
m U 

This equation is valid for both flow regions i.e. for where the pressure gradient is positive
and for where the pressure gradient is negative. The pressure distribution is obtained by
integration of Eq. 2.49, with the result

 2n  1 
n
2 R x '   x'
2 2 n 1
2  x'2  
H o  
p'    dx'
 
2 n 1 (2.51)
 n  1  x '2

The computation of the integral in Eq. 2.51 must be performed numerically and the
boundary condition p' x'    0 has been used.

It is simpler, computationally, to integrate Eq. 2.51 from λ (where p' 0 ) to x'


with the result

 2n  1 
n
2 R    x'
2 2 n 1
2  x'2  
H o x '
p'    dx' (2.52)
 n  
1  x'2
2 n 1

Then, the maximum pressure occurs at x'   , so we may write

17

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
 2n  1 
n
2 R  2  x'2  
H o  1  x'2 2 n 1
pmax '    dx' (2.53)
 n   

The dimensionless parameter λ can be obtained then from Eq. 2.51 taking into account
the boundary condition which is described by Eq. 2.43. Thus λο is found from the
following equation

  2  x'2
n 1
 
2
 x' 2 
0 dx' (2.54)

1  x' 
2 2 n 1

Note again here, that the notation λο is used to refer to that value of λ corresponding to the
entrance condition p' x    0 . Then it can be calculated based on the dependence to
the power–law index n as shown at Figure 2.4, where it is apparent that the effect of non–
Newtonian behavior is to increase the sheet thickness.

Figure 2.4. Thickness of the calendered sheet, in terms of λο or H/Ho, as a function of the power–law index
n [1].

Once λο is determined from n, the maximum pressure may be found by the


following equation

p'max  2n  1 
n
o  2
 x' 2  dx'
    1  x' 
o
(2.55)
2 2 n 1
2R H o  n   o

18

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Chapter 3

Literature Survey of Calendering


A description of the calendering process as well as the hydrodynamic behavior of
the fluid has been extensively studied and can be found in the literature. The earliest
reported study concerning the calendering process was performed by Gaskell [5]. The
analysis was based on the assumption that the diameter of the rolls is large enough
compared to the gap between the calenders, thus approximating the flow as one
dimensional. The developed analysis by Gaskell employed the flow of Newtonian fluids
as well as Bingham plastic fluids.

Middleman [1] and McKelvey [6] summarized previous works on calendering and
described a model for the calendering of Newtonian and power–law fluids. In his model
they treated the flow with the Lubrication Approximation Theory simplifying the Navier–
Stokes equations. By this, the analysis was restricted in unidirectional flow ignoring any
fluid flow in the cross–machine direction (this direction is the same with the direction of
the symmetry axis of the calenders). Although, such an approximation seems rational and
holds for processes where the width of the sheet is very large, it is far from reality when
the sheet is narrow. A theoretical approach on studying the calendering of power–law
fluids was also showed by Brazinsky et al [7].

Alston et al [8] developed a one dimensional model to study the calendering


process of non–Newtonian fluids which are described by a hyperbolic tangent model,
they called it tanh model, exhibiting Newtonian behavior at high and low shear rates and
a shear–thinning behavior at several intermediate shear rates. The authors also suggested

19

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
that extension of their work is the solution of a two–dimensional flow between the
calenders.

Kiparissides and Vlachopoulos [9] studied symmetric and asymmetric calendering


of Newtonian and non–Newtonian fluids numerically. In their model they simplified the
Navier–Stokes equations by the Lubrication Approximation Theory and the resulting
equations were solved by the finite element method. The pressure distribution along the
centerline for a Newtonian fluid showed very good agreement at lower values of λο while
at higher values there was a substantial deviation due to the coarse mesh they used
consisting of long and narrow elements. The power–law viscosity model they also used in
their analysis, produced results with the maximum pressure being substantially higher for
Newtonian fluids than the non–Newtonian ones and the asymmetry (in calenders radii
and speed) did not influence the results for the Newtonian fluid, but reduces the
maximum pressures of non–Newtonian fluids as shown at Figure 3.1. Moreover, the
same authors [10], used a finite–differences non–isothermal model to study the
temperature distribution in the calendering process due to viscous dissipation phenomena,
reporting two maxima of the temperature near the surfaces of the calenders caused by the
large amount of shear.

Figure 3.1. Pressure profiles along the centerline for calendering of a Power – Law fluid with n=0.25. The
calenders (R=100 mm, Ho=0.1 mm) are rotating at different velocities [9].

20

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Furthermore, Agassant [11,12], studied the flow between the calenders
theoretically using a two–dimensional Finite Elements method, for rigid PVC, and also
carried out experiments concerning the formed melt bank. In the computational method
used, the streamlines in the calendering melt bank where calculated and the pressure
distribution, in good agreement with the experiment, using the Lubrication
Approximation. The numerical results, concerning the streamlines, are shown at Figure
3.2 in which there are two recirculation patterns. At Figure 3.3 is shown the flow pattern
in the calender melt bank as observed experimentally with a huge vortex covering a large
portion in the melt bank, while two smaller ones develop near the entering sheet and the
nip region. It was also pointed out that when the calendering process is unstable (too
large melt bank or too high output rate) the bank as well as the recirculating flow patterns
flow in the third direction, and a spiral flow occurs. The spiral flow in the third direction
and the existence of three vortices was firtsly observed by Unkrüer [13] and more results
by his study will appear further on in the present thesis.

Figure 3.2. Numerically calculated streamlines in the melt bank [11, 12].

Figure 3.3. Experimental flow pattern in the melt bank of calendered PVC [11, 12].

21

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Vlachopoulos and Hrymak [14] studied on a theoretical and experimental
framework calendering of rigid PVC with an isothermal power–law model. The
Lubrication Approximation method of solution with no–slip gave significantly higher
pressure distribution in the nip region than the experiment. However, the experimental
pressure was measured 80 mm from the center without reporting which was the
experimental maximum pressure at the center. The Lubrication Approximation with slip
predicted lower maximum pressure than the analysis with the no–slip approximation.

Moreover, Mitsoulis et al. [15,16] performed a two–dimensional, non–isothermal


analysis with a purely viscous model without using the Lubrication Approximation
Theory for Newtonian and power–law fluids (rigid PVC with data from Vlachopoulos et
al [14]) by means of the Finite Elements method. In the analysis they captured the two
recirculation flow patterns as it can be seen at Figure 3.4 and the temperature distribution
for the region between the calenders including the melt bank as shown at Figure 3.5.
Lubrication Approximation cannot predict the asymmetric recirculation patterns in the
melt bank. The analysis also predicts the pressure distribution along the centerline. The
numerical scheme used gave reasonable predictions for the pressure distribution
especially when slip is included.

Figure 3.4. Calculated streamlines for the calendering of rigid PVC with data obtained from Vlachopoulos
and Hrymak [14], indicating two recirculation regions [15,16].

22

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 3.5. Calculated isotherms for the calendering process of rigid PVC [15,16].

More recent studies were carried out by Sofou and Mitsoulis [17, 18] for the
calendering of viscoplastic and pseudoplastic materials using the Lubrication
Approximation with no–slip and slip at the calenders surfaces but the results were of little
importance.

Moreover, concerning the processing of carbon nanotubes with calendering a


couple of studies are presented obtained from the literature [19, 20]. Kovacs et al. [19]
studied the conductivity of Multi–Wall Carbon Nanotubes (MWCNT) by the help of a
dissolver disk or additionally with a three roll calender. Specifically it was concluded that
nanotube agglomerates were efficiently separated and dispersed through the calendering
process due to the high shear forces they undergo in the nip region of the calender but
this is not favourable for electrical conductivity applications. It was also pointed out that
the flow forces exerted on the nanotubes caused attrition with a reduction of their aspect
ratio from 1000 to 60. Moreover, Seyhan et al. [20] made a study which exhibits that a
three–roll–milling process (a process which is very similar to calendering) is sufficient
for the dispersion of carbon nanotubes in a thermoset polyester resin than other processes
such as direct mixing.

In all the studies presented above, the calendering process was studied using two–
dimensional numerical schemes (mostly Finite Elements Method) or one–dimensional

23

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
models such as the Lubrication Approximation. Of course, such approximations, give
reasonable results concerning the pressure distribution along the centreline using no – slip
or even slip assumptions at the calender surfaces. However, these studies neglect any 3D
effects that take place, such as the spreading of the sheet from the entrance to the exit or
the pressure distribution in the cross – machine direction as observed by Unkrüer [13].
The only known 3D study using numerical methods is by Luther and Mewes [21] who
studied 3D effects in the melt bank when calendering non–Newtonian fluids. An
extensive experimental analysis, as stated above, was carried out by Unkrüer [13]. More
on this matter in calendering will be discussed further on in the present thesis.

24

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Chapter 4

Finite Volume Method


Introduction

For all of the calculations at the present study the OpenFOAM (Open Source Field
Operation and Manipulation) package is used which employs the Finite Volume Method
(FVM) which subdivides the flow domain into a finite number of smaller control
volumes. The transport equations are then integrated over each of these control volumes
by approximating the variation of flow properties between mesh points with differencing
schemes. The solver used in the present work is simpleFoam which is a steady – state
incompressible solver that uses the FVM to simulate the flow of Newtonian and non –
Newtonian fluids with the SIMPLE algorithm. The SIMPLE acronym stands for Semi –
Implicit Method for Pressure Linked Equations and it was adopted by Patankar and
Spalding [22] and it is actually a guess – and – correct procedure for the calculation of
pressure. This section gives only a brief overview on the FVM [23,24], for a complete
discussion on the FVM the reader is referred to An Introduction to Computational Fluid
Dynamics by Versteeg and Malalasekera [25].

4.1 Discretisation of Continuous Operations

A computer can only process discrete quantities and therefore continuous


quantities like space and time need to be discretised before processed. There are two
main distinctive discretisation procedures that must be made in order to obtain the
solution of a flow problem. These procedures are:

25

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Spatial discretisation: Discretisation of the solution domain produces a computational
mesh on which the governing equations of the flow are subsequently solved. The
procedure can be split into two parts, discretisation of space and time. The discretisation
of space for the Finite Volume Method used in this study requires a subdivision of the
flow domain into control volumes (CV) or cells which do not overlap and completely fill
the computational domain. It follows naturally that the boundary of the flow domain is
decomposed into faces (f ) connected to the cell next to it. The cell connected from the
lower side of the face is called “owner” (P) and the cell connected from the higher side of
the face is called “neighbor” (N) as shown at Figure 4.1. Generally, variables are stored at
the cell centroids although they may be stored on faces or vertices. A cell is bounded by a
set of flat faces with no limitation on the number of faces or their alignment, which can
be called arbitrarily unstructured. The capability of OpenFOAM to discretise with
arbitrary control volumes and produce unstructured meshes gives considerable freedom
in mesh generation and appears to be useful in meshing complex geometrical
configurations in three spatial dimensions.

Figure 4.1. Arbitrary control volume [24]. (P) refers to the “owner” cell, (N) refers to the “neighbor” cell,
(f) denotes the shaded area which is the face and (S) is the face area vector.

Transport equations discretisation: The partial differential equations that are


appropriate for certain cases must be discretised in order to be processed. OpenFOAM
solves certain fields, like pressure and velocity, through a translation from the partial
differential equations and the computational flow domain into a set of algebraic
expressions. Specifically, the simpleFoam solver implements the Navier–Stokes

26

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
equations and supports the option to choose between turbulent or laminar flow. Since
polymer melts are characterized by very low Reynolds numbers (i.e. 10-2 – 10-4) their
flow behavior is characterized as creeping and laminar. This was presented in a previous
chapter. For that reason, it was chosen the option of laminar flow for the solver.

4.2 Governing and Constitutive Equations of the Solver

Fluid flow, in general, is mathematically described by the conservation of mass,


momentum and energy. The simpleFoam solver, solves the equations presented below.
These equations are the conservation of mass in a fully three dimensional form and the
fully three–dimensional Navier–Stokes equation along with a constitutive equation that
relates the stresses with the rates of strain. The study of any thermal effects and
temperature distribution that take place in the calendering process is out of the scope of
the present thesis, thus the energy equation was no solved. The general form of a
conservation equation for a flow quantity  if the flow is in steady state, is:

  ( U )    (  )  S (4.1)

where ρ is the density, U= U(ux,uy,uz) is the velocity,  is the diffusivity and S is a

source term. The transport equation for the conservation of mass, or continuity equation,
is derived by setting   1 in Eq. 4.1 and not having any mass source terms. This
subsequently leads to the following equation if incompressibility is taken into account

U  0 (4.2)

The fully three – dimensional Navier – Stokes equation, neglecting any


gravitational effects, and accepting that the flow is steady in time is the following

  ( UU )  p    (4.3)

27

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
where  is the stress tensor, a second order tensor, which was showed in a previous
chapter. The equation above, was presented in a previous chapter but in a slightly
different form. In the solver the inertial term of the momentum equation is taken into
account. In the following sections is presented how each term of the Navier–Stokes
equation is discretised.

4.3 Discretised Equations

4.3.1 System of Equations

The discretisation and linearization procedure outlined in the previous sections


produces a linear algebraic equation for each control volume. The exact form of these
linear algebraic equations depends on the governing equation and the discretisation
practices used, but they can be re–written in a generic formulation, in matrix form such as
below

   b (4.4)

where [A] is a tridiagonal matrix in the form below

  C 
 C  C 
 
 C  C  (4.5)
 
C  C
 
 C  C
  C  

[x] is the column vector of dependent variable, in the form below

 1 
 
  2 (4.6)

 
  

28

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
and [b] is the source vector

 b1 
b 
b 2 (4.7)

 
b 

The solution of this system approximates the solution to the original equations at
some pre–determined locations in space (and time for unsteady flows). The numerical
methods for the solution of the above type of matrix fall into two main categories, direct
methods and iterative methods. Direct methods involve direct solution of Eq. 4.4 and
require inversion of the matrix [A]. Iterative methods, start with an initial guess and then
continue to improve the current approximation of the solution until some tolerance is met.
For direct methods, the number of operations necessary to reach for the desired solution,
scales approximately with the cube of the number of equations (or the unknowns), thus
making them computationally very expensive especially for large systems. On the other
hand, the total number of operations required by iterative methods, typically of the order
of N per iteration cycle, cannot be predicted in advance. Stronger still, it is not possible to
guarantee convergence unless the system of equations satisfies fairly exacting criteria.
The main advantage of the iterative methods is that only non–zero coefficients of the
equations need to be stored in memory thus, making these methods for solution more
economical computationally.

On a specific number of iterations, for iterative methods, the solution will be  n ,


and it will not fully satisfy Eq. 4.4, in the sense of a non–zero residual rn.
Mathematically, this can be expressed as

 n  b  r n (4.8)

The aim is to minimize the value of the residual. If that value is small enough, then it can
be approximated that

29

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
 n    (4.9)

In the present study the conjugate gradient method is used in order to reach for the
solution after a certain number of iterations which was introduced by Hestenes and
Stiefel [26].

4.3.2 Linearization of the Navier – Stokes Equations

In the present study each term of the Navier – Stokes equation (Eq. 4.3) is
discretised. Finite volume (FV) discretisation of each term is formulated by integrating
the terms of the equations over a cell volume V (control volume). If not all, at least most,
spatial derivatives are converted to integrals over the cell surface S that bounds the
volume by means of Gauss’s theorem

    dS
V s (4.10)

where S is the surface area vector and  can represent any variable. The calculated
values of pressure (p) and velocity (U) are stored at the same nodal points which are
actually the centroid of each cell. Such an arrangement of the computational domain is
called collocated arrangement. OpenFOAM utilizes such type of arrangements, like in
the present study.

The Diffusive Operator


The diffusion term is integrated over a control volume and linearised as follows

   (U)dV   UdS    f S f  U  f (4.11)


V S
f

where μ is the viscosity for a Newtonian fluid and the index f refers to the face
interpolated values. The previous can be discretised when the length vector d between the

30

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
centre of the cell of interest P and the centre of a neighboring cell N is orthogonal to the
face

U  UP
S f  U  f  S f (4.12)
d

In the present work, the term was evaluated using Gaussian integration which requires
the interpolation of values from the cell centroid to the face centers. The interpolation
scheme used is central differencing for the μf and the surface normal gradient U , which
is second order accurate. As explained in detail in Breuer [27], a second–order accurate
scheme is even appropriate for LES of turbulent flows and consequently is reliable for the
simulation of laminar flows carried out within the present study. The same scheme was
also used by Biswas et al. [28] who investigated laminar and creeping flows in a
backward–facing step at low and moderate Reynolds numbers.

The Convective Operator


The convection term, as previously, is integrated over a control volume and linearised as
shown below using Gaussian integration

   ( UU )dV   dS  UU    S  U


V S
f
f f U f  FU f
f
(4.13)

where F in Eq. 4.13 represents the mass flux through the face. The face field U f can be

evaluated using a variety of schemes. In the present work, the used interpolation scheme
is the upwind differencing which is first order accurate but bounded and determines  f

from the direction of flow. That scheme was used for stability reasons.

The Gradient Operator


The gradient term, which is described here, is an explicit term. Usually it is evaluated
using Gaussian integration (as in the present study) which requires the interpolation of
values from the control volume centroid to the face centers. The interpolation scheme

31

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
used is central differencing. Therefore, applying Eq. 4.10 to the volume integral we have
the following

 pdV   dS  p   S
V S
f
f pf (4.14)

where p f is the value obtained from the interpolation and Sf is the area of the surface

that bounds the control volume. There are more ways to evaluate the gradient term. For
more information the reader is referred to the OpenFOAM Programmer’s Guide [29] and
the User’s Guide [23].

4.4 Numerical Setup

For the simulations, the simpleFoam solver of the OpenFOAM package was
employed which is a steady – state solver for incompressible, viscous fluids. The flow is
assumed to be laminar. For the pressure calculation, a pressure–correction equation
taking mass conservation expressed in Eq. 4.2 into account was solved iteratively along
with Eq. 4.3. The entire procedure follows the well–known SIMPLE algorithm proposed
by Patankar and Spalding [22]. The solution of the algorithm was met until a desired
tolerance was reached. It was assumed the tolerance for the pressure and the momentum
equation to be 10-7 both for the Newtonian and the Generalized Newtonian Fluid.

32

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Chapter 5

Calendering Validation with OpenFOAM

Introduction
In order to validate the simpleFoam solver against the approximated analytical
solutions, a computational flow domain is constructed based on an arbitrary geometry
which is rather common in the commercial calendering machines. Firstly, it was tested a
geometry with a Newtonian fluid and further on with the Generalized Newtonian Fluid
model for the same geometry. The approximate geometry along with the boundary
conditions of the Lubrication Approximation Theory, pose a two–dimensional solution
but on a three–dimensional geometry. For the sake of saving computational time, it could
be constructed a two–dimensional geometry. However, the present three–dimensional
geometry will be, in subsequent chapter, the basic computational flow domain on which
the three dimensional effects of calendering will be studied after small modifications in
the flow domain.

5.1 Newtonian Fluid

5.1.1 Flow Domain Geometry

At this chapter the flow behavior of a simulated Newtonian isothermal fluid


between the calenders is studied using OpenFOAM and is directly compared to the
analytical solution. The flow domain generated with the software is presented at Figure
5.1 and Figure 5.2 and the geometry dimensions are shown at Table 5.1. The x–axis will

33

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
represent the machine direction (that axis will always be parallel to the centerline shown
at Figure 2.1). The y–axis will represent the vertical distance from the top to the botton
surface of the calenders and finally z–axis will represent the cross – machine direction
(which means that z–axis is parallel to the symmetry axis of each calender). The values
for the geometry used are arbitrary but close enough to the dimensions of commercial
calenders.

Figure 5.1. Side view (cross–section) of the computational flow domain, for xi,res. = -100.00 mm

Figure 5.2. Inclined view of the computational flow domain for xi,res. = -100.00 mm.

34

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
The results of the present simulations are based on the pressure distribution along the
centerline (x–direction), which is the machine direction, and then they are directly
compared to the analytical solution from the Lubrication Approximation Theory. The
analytical solution was obtained using the CALENDERCAD [29] which solves the
Lubrication Approximation equations with a fourth–order Runge–Kutta formula and can
calculate the pressure distribution along the centerline as well as the approximated point
where the sheet detaches from the calenders with no–slip.

Table 5.1.
Dimension of the Geometry Value
Nip gap (2Ho) 1.000 mm
Final sheet thickness (2Hf) 1.226 mm
Roller diameter (D) 400.000 mm
Roller velocity (V) 0.0334 m/s
Roller width (W) 500.000 mm
Attachment point from -40.000 mm
Lubrication Approx. (xi)
Detachment point (xf) 6.72 mm
Infinite reservoir point (xi,res.) -100.000 mm
Flow rate (Q) 55 kg/hr
Viscosity (μ) 15000 Pa·s
Density (ρ) 1000 kg/m3

As far as the detachment point is concerned, the adopted value was the one
derived from Eq. 2.42 of the Lubrication Approximation theory with λo = 0.475 being the
dimensionless detachment point (or xf =6.72 mm) and Hf / Ho = 1.226 as it seems a very
reasonable approximation. The geometry of the computational flow domain is calculated
through the following simple equation by which the shape of the calenders is cylindrical
(and thus, their cross–section circular)


 x   H o  R  R 2  x 2 
1
2
(5.1)

where H(x) is the half gap between the rollers at each position x. However, Middleman
[1] refers that for positions which are relatively close to the nip region, in order to

35

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
simplify his analysis with the Lubrication Approximation, the surfaces of the calenders
can be approximated as a parabolic shape. Thus Eq. 5.1 is approximated as

 x2 
x   H o 1 
  (5.2)
 2 H o 
R

Since, in the present study an infinite reservoir is assumed upstream which is located far
behind the nip region, the only reasonable approximation for the construction of the
computational domain was to treat the calenders surface as circular and use 5.1. A similar
procedure is followed by Mitsoulis [17,18] for the numerical simulation of calendering
viscoplastic fluids using the Finite Element Method. Following Middleman’s [1]
approach it is also assumed in the present study that far behind the nip region the feeding
is done from an infinite reservoir. It was also approximated that pressure becomes zero
(p=0, x→-∞) far away behind the nip region. The computational domain shown at
Figures 5.1 and 5.2 does not necessarily represent the true physical boundaries of the
process, since in the real process the feeding is done with a sheet of finite thickness and a
melt bank appears near the attachment point. The shape of the free surface that describes
the melt bank is however, beyond the scope of this study. The approach of the
computational domain adopted in the present study is done for the sake of imposing the
boundary conditions as described by the Lubrication Approximation and for further
comparison with the analytical solution.

5.1.2 Boundary Conditions

The flow domain consists of regions that coincide with the boundaries of the
physical domain. To correctly simulate the conditions between the two rollers, boundary
conditions at all boundaries also need to be specified. Below the used boundary
conditions of the Newtonian case study are shown following the same notation as in
[23,30] and the boundaries are shown at Figure 5.3. In the formulation below the vector n
denotes the direction normal to the boundary.

36

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
 Top and bottom rollers
p
- Pressure: Normal gradient of pressure is zero, 0
n
- Velocity: no – slip at the roll surface (tangential velocity ut=V, normal
velocity un=0)

 Sides
p
- Pressure: Normal gradient of pressure is zero, 0
n
U
- Velocity: Normal gradient is zero, 0
n

 Outlet
- Pressure: Fixed value, p = 0
U
- Normal gradient of velocity is zero, 0
n

 Inlet
- Pressure: Fixed value, p = 0
- Velocity: Known flow rate Q    ut  W  2Hi 
where 2Hi is the distance between the rollers at the inlet.

Figure 5.3. Side view of the computational domain with notation of each boundary.

37

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Some comments concerning the boundary conditions should be made to this
point. First of all it was approximated that the inlet boundary is far away from the nip
region of the calenders and the material is fed continuously from an extruder or any other
device. The fluid flow at the inlet is not delayed by any means of obstacles (i.e. there is
not any pressure build–up) thus, the only reasonable assumption that can be made is to set
the pressure equal to zero as this boundary is in contact with the atmosphere. Concerning
the velocity boundary condition at the same boundary, it is assumed that the flow rate is
known with a constant velocity profile. For the velocity boundary at the exit, where the
detachment point lies, it is accepted that the gradient of the velocity normal to the
boundary is zero, as presented above. This implies fully developed flow along the normal
direction. On that same boundary, it was imposed that p=0. In Lubrication
Approximation according to numerous publications such as Middleman [1], McKelvey
[6], Mitsoulis [15-18] and Tanner [31] the boundary condition imposed at the detachment
point is p=dp/dx=0. In a similar, almost convenient manner, and since that boundary is
actually a free surface, it is treated as a surface with p=0 (as it is assumed that is in
contact with the atmosphere) and zero gradient for the velocity, considering that in the
present study the detachment point is obtained from the Lubrication Approximation. By
this it is implied that from this boundary there is outflow. Concerning the side
boundaries, it was approximated that there is no spreading of the sheet in the z–direction
(cross–machine direction) and the prescribed boundary conditions basically assume
calenders with infinite length along that direction and in that the effects along the z–
direction are not taken into account for the present chapter.

5.1.3 Results and Discussion

The results obtained from simulations for the pressure distribution along the
centerline are presented at Figure 5.4 for a Newtonian fluid with the well know bell–
shaped curve. As it was previously discussed, it was assumed that far behind the nip
region there is an infinite reservoir of the material where it is assumed that p=0 as
described by the Lubrication Approximation Theory. The question that arises here is how
much far that point should be. At first it was assumed that p=0 at xi =-40.00mm (value
obtained from CALENDERCAD [29]) and at the second simulation that p=0 at xi,res.=-

38

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
100.00 mm as shown at Figure 5.4. Also, based on the analytical solution by
CALENDERCAD [29] at the attachment point the pressure is not taken to be zero but it
has a small pressure build–up. This might be due to the approximations the software
performs to calculate the attachment point. So, the location of xi,res.=-100.00mm seems as
a rational approximation as the numerical results exhibit very good matching with the
analytical solution. The maximum pressure reported by the Lubrication Theory is
Pmax=10.77MPa and the simulation gives Pmax=10.75MPa. The simulations were
performed for a flow domain filled with 225000 cells and the detachment point was at
6.72 mm and calculated from the Lubrication Approximation Theory.

Lubrication Approximation
OpenFOAM, xi,res.=-40.00 mm
12
OpenFOAM, xi,res.=-100.00 mm

10

8
Pressure (MPa)

0
-0,11 -0,10 -0,09 -0,08 -0,07 -0,06 -0,05 -0,04 -0,03 -0,02 -0,01 0,00 0,01 0,02
Position (m)

Figure 5.4. Pressure distribution along the centerline with p=0 at xi=-40.00 mm and at xi,res.=-100.00 mm.
The pressure is plotted at W/2 (middle section).

Furthermore, the behavior of the numerical solution with different levels of


refinement of the computational flow domain was investigated as shown at Figure 5.5.
For 2000 cells the solution shows a relatively small oscillation (wiggles) around the point
of the maximum pressure and this is caused due to the poor domain discretisation,
whereas for the 27750 cells the numerical solution was good enough to the one of 225000
cells.

39

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Lubrication Approximation
12
OpenFOAM, Cells=225000
OpenFOAM, Cells=27750
OpenFOAM, Cells=2000
10

Pressure (MPa)
6

0
-0,10 -0,09 -0,08 -0,07 -0,06 -0,05 -0,04 -0,03 -0,02 -0,01 0,00 0,01 0,02
Position (m)

Figure 5.5. Pressure distribution along the centerline for three different levels of mesh refinement at W/2.

At Figure 5.6 the effect of displacing the detachment point 0.2 mm and 1 mm is
exhibited. In the curves of Figure 5.6 we show how the displacement of this point effects
the shape of the pressure distribution along the centerline as well as its peak value. This is
done since in practice, even Newtonian fluids tend to slip at the exit of a die or a channel
[34,35], which means that the exact point or line of detachment is difficult to predict. The
analytical solutions derived in a previous chapter, do not take into account any slip
phenomena that occur at the surfaces of the calenders and the value for the detachment
point is obtained through various simplifications which were presented earlier in the
present section. At first, we displace the detachment point by approximately 0.2 mm
backwards and then at 0.2 mm forward. The shape of the bell–shaped curve is not altered
at all with the maximum pressure preserved to Pmax=10.75 MPa approximately. By this, it
can be concluded that the exact location of the detachment point is not a critical
parameter in estimating the maximum pressure as well as the shape of the pressure curve.
For larger displacements of the detachment point, of the order of 1 mm backwards and
forwards, the plotted pressure profile had the same distribution as in the previous
simulations and the maximum pressure was Pmax=10.75 MPa for the displacement at xf
=5.72 mm and Pmax=10.71 MPa for the displacement at xf =7.72 mm. Therefore, the above

40

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
conclusion can be extended for uncertainty in the location of the detachment point of the
order of ±1 mm.

Lubrication Approximation
OpenFOAM, xf = 6.72 mm
OpenFOAM, xf = 6.52 mm
OpenFOAM, xf = 6.92 mm
12,0
OpenFOAM, xf = 5.72 mm
10,5 OpenFOAM, xf = 7.72 mm

9,0

7,5
Pressure (MPa)

6,0

4,5

3,0

1,5

0,0

-1,5

-0,10 -0,08 -0,06 -0,04 -0,02 0,00 0,02


Position (m)

Figure 5.6. Pressure distribution concerning the displacement of the detachment point. The Lubrication
Approximation gives a detachment point at xf =6.72 mm. The infinite reservoir is assumed to start at xi,res.=
-100.00 mm. All the values are plotted at W/2 (middle section).

At Figure 5.7 and Figure 5.8 the velocity field and the streamlines are shown
respectively for a flow domain of 760000 cells and for xi,res.=-100.00 mm and xf =6.72
mm. The numbers of cells for that simulation was increased to 760000 cells because we
were interested to capture the details of the flow regime. The results show clearly two
symmetric to the centerline recirculation patterns, something that has been mentioned in
numerous publications discussed in a previous chapter for calendering of Newtonian
fluids. These vortices are a result of the drag induced flow of the Newtonian fluid
between the calenders and the back–flow due to the pressure build–up at the nip region.
From Figure 5.9, it can be concluded that the maximum velocity of the flow is developed
in the nip region. Due to the fact that the distance between the calenders in the nip is the
minimum, the velocity of the material increases. Such velocity behavior is observed in
entry flows where the material enters to a small tube from a large reservoir. Extensive
simulations on the entry flows were made by Mitsoulis [32] in his doctorate thesis.

41

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 5.7. Streamlines for a flow domain of 760000 cells. xi,res. = -100.00 mm and xf =6.72 mm.

Figure 5.8. Variation of fluid speed along the flow domain of 760000 cells. xi,res. = -100.00 mm and xf =
6.72 mm.

Furthermore, it can also be concluded that the outlet boundary (which is marked in Figure
5.3) appears to have the same velocity with the velocity of the calenders implying that the
Newtonian fluid leaves the point of detachment with the same velocity.

42

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 5.9 Velocity vectors for a flow domain of 760000 cells. xi,res. = -100.00 mm and xf = 6.72 mm. The
vectors are plotted at a computational domain slice located at W/2.

Middleman [1] reports such a behavior in his study with the analytical approximation,
where it is assumed that after the detachment point the velocity is constant (plug flow)
and equal to the velocity of the calenders. This is validated and further on the present
study by calculating the velocity profile at the exit boundary. However, in Mitsoulis [33]
study it is reported that the material at the exit leaves the rollers with a little lower
velocity. Finally, the velocity vectors are plotted on a slice of the computational domain
in the middle region (located at W/2) as shown at Figure 5.9.

Moreover, at Figures 5.11, 5.12 and 5.13 the velocity components and their
magnitude in the y–direction are plotted for different positions as showed at Figure 5.10.
It can be seen at Figure 5.11 that at x=-70.00 mm, ux is maximum near the calenders and
near the centerline decreases, taking negative values. The negative values imply backflow
and this is why the recirculation patterns occur. The velocity component in the y–
direction, uy, appears to be maximum near the surfaces of the calenders and holds small
to zero values for a region close enough to the centerline. This implies, that near the
centerline the dominant velocity component is in the x–direction. The velocity in the z–
direction, uz, is approximately zero something which was anticipated based on the set of
boundary conditions at the two side boundaries. The velocity magnitude exhibits
maximum values near the surfaces of the calenders, and as we move to the centerline

43

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
from the calender surfaces, it decreases and finally exactly on the centerline increases
again. The increase on the centerline is due to the occurring backflow and since the
velocity magnitude is always a positive value, its values lie on the positive side of the
axis, resulting in Umag. profile showed at Figure 5.11. Similar velocity profiles are found
at x=-30.00 mm, with the maximum absolute value of ux being smaller than the location
of x=-70.00 mm.

Figure 5.10. Cross–sectional area with the velocity magnitude field and velocity vectors with indications at
four different arbitrary locations normal to the xz plane.

ux
uy
0,215
uz
0,210 Umag.
y - direction (m)

0,205

0,200

0,195

0,190

0,185
-0,02 -0,01 0,00 0,01 0,02 0,03 0,04

Velocity (m/s)

Figure 5.11. Velocity components profiles and velocity magnitude profile at x=-70.00 mm.

44

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
This means that as the fluid progresses to the exit the recirculation effects are becoming
less intense, something which can be seen at Figures 5.13 and 5.14. At Figure 5.13 the
profiles of ux and Umag. match, thus the dominant flow occurs in the x–direction and since
the velocity profile is parabolic with positive values, there is not any recirculatory
patterns at that location. The velocity components at the other directions are
approximately zero. For the exit boundary, at xf =6.72 mm the velocity profile exhibits
small variations but it can be considered as flat as it can be seen at Figure 5.14. This is
the result of the velocity applied boundary condition at the exit, by which the velocity
gradient normal to the boundary is zero and as in the previous location ux and Umag.
curves match. The matching of the curves implies that at the nip region as well as after
the nip the flow can be considered as one dimensional as the dominant velocity
component lies in the x–direction. The present results, validate the behavior of the
velocity presented by previous studies, namely in Middleman [1], Gaskell [5] and
Mckelvey [6].

ux
0,204
uy
0,203 uz
Umag.
0,202
y - direction (m)

0,201

0,200

0,199

0,198

0,197
-0,010 -0,005 0,000 0,005 0,010 0,015 0,020 0,025 0,030 0,035
Velocity (m/s)

Figure 5.12. Velocity components profiles and velocity magnitude profile at x=-30.00 mm.

45

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
ux
0,2012
uy
0,2010 uz
Umag.
0,2008

y - direction (m)
0,2006

0,2004

0,2002

0,2000

0,1998
0,00 0,01 0,02 0,03 0,04 0,05
Velocity (m/s)

Figure 5.13. Velocity components profiles and velocity magnitude profile at x=0 (nip region).

ux
uy
0,2012
uz
0,2010
Umag.
0,2008
y - direction (m)

0,2006

0,2004

0,2002

0,2000

0,1998

0,000 0,005 0,010 0,015 0,020 0,025 0,030 0,035


Velocity (m/s)

Figure 5.14. Velocity components profiles and velocity magnitude at x=6.72 mm (exit boundary).

46

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
5.2 Generalized Newtonian Fluid Model

5.2.1 Flow Domain Geometry and Boundary Conditions

At the present chapter the pressure distribution along the centerline for a non–
Newtonian fluid is investigated neglecting any thermal phenomena that physically take
place in the process. Thus the energy equation is not solved in the present study. The
validation of the simulations is done based on the work of Vlachopoulos and Hrymak
[14]. For that reason, and again assuming that there is an infinite reservoir far behind the
nip region, we set the point of xi,res.=-100.00 mm which will be the location of the inlet
boundary for the computational flow domain. In the previous Newtonian study the
equations of the Lubrication Approximation Theory are solved analytically and the value
of the detachment point is calculated exactly. However, for the power–law fluid the
Lubrication Approximation is applied through numerous approximations and the value of
the detachment point is more difficult to be calculated analytically as with the power–law
viscosity model an extra parameter is inserted in the equations (that is power–law index
n). A recent publication of Mitsoulis [36] concludes that the exact location of detachment
(and attachment point) in calendering still remains an unresolved issue. The detachment
point according to the Lubrication Approximation is obtained analytically to be at xf
=4.37 mm with the assumption that the feeding is done from an infinite reservoir far
upstream. In the present study, we accept xf =4.37 mm as a sufficient initial guess of the
detachment point. The feeding of the material is approximated to be through an infinite
reservoir that comes off an extruder or any other continuous polymer processing
machinery. Concerning the viscosity model, the power–law viscosity is utilized (Eq.
2.23) with m=46120 Pa·s and n=0.34 as shown at Table 5.2 below. Although there is an
exponential or Arrhenius type dependence of viscosity with the temperature, in the
present study we assume that the viscosity is no affected by temperature. Moreover, as in
the simulations for the Newtonian fluid, the surfaces of the calenders are assumed to be
cylindrical (Lubrication Approximation Theory assumes parabolic calenders surfaces)
and the polymer cannot slip. Since we examine the pressure distribution only in the x–
direction, it seemed rational to utilize the boundary conditions of the previous study with
the Newtonian fluid.

47

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
5.2.2 Results and Discussion

The first simulation performed was to test the geometry given in Table 5.2 with a
Newtonian fluid and the detachment point was known a priori from the Lubrication
Theory to be xf =4.37 mm. The maximum pressure calculated by the Lubrication
Approximation was pmax =130.15 MPa while the simulation gave a maximum pressure of
pmax.=129.5 MPa in excellent overall agreement. Again here, the boundary conditions of
the flow domain are the same with the Newtonian case described in a previous chapter.
The second simulation employed the same geometry as above but the viscosity was
described with the power–law model. The pressure distribution with the simulation is
shown at Figure 5.15. From a very first glance it can be seen that OpenFOAM over–
predicts the maximum pressure of pmax.=8.51 MPa, with respect to the analytical solution
derived by Vlachopoulos et al [14] in which pmax.=7.1 MPa.

Table 5.2

Calendering geometry from [6]


R = 125 mm Ho = 0.3 mm W = 450 mm xf = 4.37 mm
Hf = 0.74 mm Vroller1 = 0.07383 m/s Vroller2 = 0.07366 cm/s Q =100 kg/hr
Viscosity parameters from [6], for power law viscosity model
m = 46120 Pa·s n = 0.340 ρmelt = 1300 kg/m3

Since this might have been due to the poor domain discretisation, simulations
were performed in which the number of cells was increased to 360000 and 760000 and
the results are shown at Figure 5.16. Mesh refinement did not affect the behavior of the
predicted pressure distribution as it still over–predicts the maximum pressure obtained
from the analytical solution from the Lubrication Approximation. Such a behavior of the
pressure distribution might occur due to the power–law viscosity model used, which is a
two parameter model and subject to limitations when used in a 2D or 3D geometry, as
explained earlier. Contrary to the predictions of the power–law model, polymer melts and
liquids, at sufficiently low shear rates, exhibit a plateau viscosity which is significantly
lower than the viscosity predicted by the power–law model.

48

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Lubrication Approximation
OpenFOAM (130000 cells)
9

Pressure (MPa)
5

-1
-0,10 -0,08 -0,06 -0,04 -0,02 0,00 0,02
x - direction (m)

Figure 5.15. Pressure distribution along the x – direction for the computational flow domain filled with
130000 cells and analytical solution with Lubrication Approximation.

Lubrication Approximation
OpenFOAM (130000 cells)
9 OpenFOAM (350000 cells)
OpenFOAM (760000 cells)
8

6
Pressure (MPa)

-1
-0,10 -0,08 -0,06 -0,04 -0,02 0,00 0,02
x - direction (m)

Figure 5.16. Pressure distribution for different number of cells filling the computational flow domain.

This means that at low shear rates they behave as Newtonian fluids. As the shear rate
increases, in most polymeric fluids, the viscosity drops (this behavior is often called
shear thinning), and on double–logarithmic coordinates one often observes nearly
straight line behavior at high shear rates. Such a behavior is shown as an example in
Figure 5.17 for HDPE. At the same figure it can be concluded that at very low shear rates
(i.e.  1 s 1 ) the power–law model severely overestimates the viscosity (when   0

49

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
the viscosity tends to infinity (    )). For that reason, the Power – Law model may
give erroneous results if used in a flow domain where regions of very low shear rate
exist. If someone is to use the power–law model he needs to truncate the viscosity at the
very low shear rates region to capture the Newtonian behavior of the material. Obviously,
there is no advantage in choosing the simple power–law model whenever finite elements,
finite volumes or finite differences are involved and for that reason it should be avoided.

Figure 5.17. Viscosity versus shear rate for High Density PolyEthylene (HDPE) at different temperatures.
The points are measured values, the solid lines represent the fit with the Carreau model and the dashed ones
with the Power – Law model [2].

Lub. Approx. with no-slip


Lub. Approx. with slip (Vlachopoulos et al )
OpenFOAM with Carreau
9 Experiment (Vlachopoulos et al )

6
Pressure (MPa)

-1
-0,10 -0,08 -0,06 -0,04 -0,02 0,00 0,02
x - direction (m)

Figure 5.18. Pressure distribution along the centerline with the Carreau viscosity model. Results obtained
with OpenFOAM, Lubrication Approximation with no – slip, Lubrication Approximation with slip [14] and
experimental measured values [14].

50

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
To predict correctly the pressure distribution between the calenders the Carreau viscosity
model is employed which is a four parameter model and can predict the viscosity at low
shear rates as it can be seen from Figure 5.17. Vlachopoulos and Hrymak [14] in their
work do not report any measured viscosity versus shear rate plots nor where the viscosity
truncation for the power–law model was carried–out. In the present study it is
approximated that the truncation of the viscosity was at 46120 Pa·s, which means that
below 1 s-1 of shear rate, the fluid behaves as Newtonian. Thus, in order to employ the
Carreau model we approximate that η0=mpower-law=46120 Pa·s, for the power–law index
n=0.34 and λ=1.25 s. With such an approximation the Carreau model can fit very well the
power–law at the high shear rate region near the nip region. The pressure distribution in
the machine direction (centerline) with the Carreau viscosity model is shown at Figure
5.17. The maximum pressure reported by the Lubrication Approximation is 7.1 MPa and
the maximum pressure OpenFOAM calculates is 7.4 MPa which is in general in good
agreement. The maximum pressure reported by the experiment is 4.02 MPa. However,
the value of 4.02 MPa is measured by a pressure transducer which was located 80 mm
from the center of the calenders as Vlachopoulos et al report and this is a reason why
such a deviation exists. Moreover, there was no report in their work for the experimental
pressure distribution for the central region along the centerline. The other reason is that
the present study has not taken into account any slip phenomena which may take place at
the surface of the calenders. The slip velocity at the surfaces of the calenders can be
measured. If slip at the calenders is inserted in the Lubrication Approximation then the
maximum pressure drops to the value of 5.60 MPa as Vlachopoulos et al [14] report in
their work. However, in the present study slip phenomena where not taken into account.

The numerical experiments carried out for the Newtonian fluid showed excellent
agreement with the analytical solution from which the detachment point was calculated
and accepted as a sufficient initial guess. However, on an effort to validate the numerical
and the experimental work by Vlachopoulos et al. [14], our results exhibited a small
deviation from the predicted and the experimental values. This was probably because
Vlachopoulos et al. [14] do not report where the viscosity truncation was done or any
measured viscosity versus shear rate curves. Since, the power–law viscosity model

51

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
exhibits a peculiarity at very low shear rates, the Carreau viscosity model was the next
rational choice, but not the only as someone could use the Cross viscosity model.
Provided that we had available viscosity versus shear rate data, the results could be
refined and correctly predict the pressure distribution.

52

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Chapter 6

3D Effects in Calendering
Introduction

The problem of rolling of sheets, such as metals and ceramic pastes is an


important problem and has been studied extensively using elastic–plastic simulation
methods to study two- and three- dimensional flows [37–45]. An interesting approach to
study spreading phenomena has been carried out by Matsumiya and Flemmings [46],
using the Lubrication Approximation Theory for the production of a metal alloy strip
from its molten state. These authors described the rheological behavior of the alloy as a
power–law fluid and assumed that the pressure did not vary in the machine direction and
the only present velocity is the one along the cross–machine direction for the region
between the rollers. Thus the pressure varies only in the cross–machine direction giving a
quasi–parabolic pressure profile. Their results are shown at Figure 6.1. Moreover, Sezek
et al. [47] in their work studied cold and hot metal plate rolling based on a three–
dimensional mathematical model considering lateral spread of the plate. They showed
that by increasing the thickness reduction of the plate the rolling force increases but
decreases with the process temperature. Three dimensional effects in the calendering of
thermoplastics have been studied by Unkrüer for the first time in his doctoral thesis [13].
He carried out experiments using a calender having two rolls of 80 mm diameter and
1600 mm length calendering PVC and PS. His conclusions were that (i) the calendered
sheet spreads and (ii) vortices are formed in and near the melt bank as shown

53

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.1. Pressure distribution in the cross–machine direction as calculated by Matsumiya and
Flemmings [46].

schematically in Figure 6.2. Also, it was observed that in the melt bank the material
conveys to the sides through a spiral motion as also shown at Figure 6.2. He also
observed, as shown in Figure 6.3, that the pressure drops to zero in the cross–machine
direction with the distance from the middle section (centerline) thus resulting in a quasi–
parabolic pressure profile as shown at Figure 6.4.

Figure 6.2. Schematic diagram of the vortices in the calendering process as observed by Unkrüer [13]. Vw1
and Vw2 are the velocities of the bottom and top roller respectively, A is the point of detachment and h A is
the final sheet thickness.

54

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.3. Spiral transport of the material, in the melt bank, from the middle region to the sides as
observed by Unkrüer [13].

Figure 6.4. Pressure distribution at various cylinder axial positions and in the cross – machine direction
near the nip region for rigid PVC as observed by Unkrüer [13]. The roll spread is 50 mm/s for both rolls,
roll temperature 185oC, minimum gap of 0.6 mm, roll diameter 300 mm and width 500 mm.

This is how Tadmor and Gogos describe Unkrüer’s observations: “The melt accumulates
in the center zone of the nip area and simultaneously undergoes flow into the nip and
sideways. The drag–induced flow leads to pressure buildup, which inevitably produces
pressure gradients in the machine and cross–machine directions”. The only known
attempt of a 3D simulation of the calendering process was carried out by Luther and
Mewes [21]. These authors however did not report any simulation of spreading or the
existence of a third vortex (two vortices are predicted by the 2D analysis as explained
earlier) or the existence of a quasi–parabolic pressure profile in the direction of the
calenders axes (cross–machine direction). Thus, their main result is the spiral motion in

55

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
the cross – machine direction as shown at Figure 6.5 which was observed experimentally
before by Unkrüer [13].

Figure 6.5. Three dimensional simulation results for the calendering of a power – law fluid. In the figure it
is shown the computational flow domain along with the spiral motion occurring in the melt bank in the z –
direction [21].

Moreover, Levine et al. [48] took into account the flow in the cross–machine
direction to predict the spreading of a power–law fluid utilizing the simplifications of the
Lubrication Approximation. They proposed a two dimensional model to describe the flow
field which seems to be similar to the Hele–Shaw approximation model for laterally
spreading flows. In the results pressure contours were calculated in the cross machine
direction, with the pressure being maximum in the central region of the sheet, but no
quasi parabolic pressure profiles were presented. It was also predicted that the sheet
spreads more when the nip region between the calenders decreases. However, in his
study, there is no reference of the work by Matsumiya and Flemmings [46].

In the present work a different approach is followed using fully 3D Navier–Stokes


equations. The spreading of liquids can be found in some publications by Scriven et al.
[49,50] and by Li et al. [51]. From a general point of view, surface tension, is an
important property of fluids and in simple words, is the force that keeps the fluid’s
molecules together. Also, the viscosity of polymer melts is roughly a million times bigger

56

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
than the viscosity of a fluid like water (μpolymers=103–105 Pa·s for polymer melts at usual
processing conditions, μwater=10-3 Pa·s). However, the surface tension of polymers is
about 1/3 of the surface tension of water (e.g. γHMWPE =29.4 mJ/m-2, γPMMA =31.7 mJ/m-2,
γwater=72.75 mJ/m-2) and for that reason, in polymer melts, surface tension effects are
rather small as compared to the viscous effects. In the real process of calendering, as the
viscous polymer melt passes through the calenders, it is squeezed. Simultaneously, and
very close to the nip region, it spreads in the cross–machine direction to a finite width.
Therefore, from a physical point of view, the macromolecular chains in polymer melts
are hold together (i.e. the material does not spread infinitely) due to a property that
primarily describes solid materials, which is called elasticity. Polymers exhibit a dual
behavior, partly as liquids (viscous) and partly as solids (elastic) as described by Ferry
[52]. In the case of calendering, the partially elastic nature of the molten polymer does
not allow it to spread infinitely. The Navier–Stokes equations do not employ any terms to
account for the elastic nature of the material since they describe the flow of liquids. For
that reason, the solution of these equations cannot predict how much the polymer will
spread in the cross–machine direction. However, if it assumed that the polymer melt fully
behaves as a fluid with a very high viscosity and neglect its elastic nature, the
applicability of the Navier–Stokes equations is valid. Then, the spreading of the molten
material is not calculated from the equations themselves but it is guessed before the
solution.

In the following section we describe the flow between the rollers with a purely
viscous model (Newtonian and GNF model) so that to investigate the behavior of the
flow at the machine and the cross–machine direction. The spreading is determined
approximately based on an initial educated guess as it will be explained further on the
present study. This assumption is done for the construction of the computational flow
domain. A few simulations, as it will be explained, can improve the initial guess of the
geometry for the flow domain. The choice of the boundary conditions for pressure and
velocity as well as approximate flow domain geometries will be discussed in the
following section with a number of simulations for a Newtonian fluid and further on for a
Generalized Newtonian Fluid employing the Carreau viscosity model.

57

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
6.1 Numerical Results for a Newtonian Fluid

6.1.1 Geometry of Computational Flow Domain and Boundary Conditions

As a first simulation it was chosen, for the geometry of the grid, the experimental
dimensions measured by Vlachopoulos and Hrymak [14]. The grid is shown at Figure 6.6
and at Figure 6.7 from another perspective. It was approximated that, for the present case,
the geometry variation (width of sheet) from the entrance to the exit is approximated to
increase linearly from the entrance to the exit, simulating the spreading of the sheet
which occurs during calendering. As described in a previous chapter the roll radius is
R=125.00 mm, the half gap at the nip region (x=0) is Ho=0.30 mm, the final sheet
thickness is 2Hf =0.74 mm, the point of detachment which was derived from the
Lubrication Approximation is xf =4.37 mm, the width of the sheet at the entry is
Wen.=225.00 mm, the width of the sheet at the very exit is Wex.=370.00 mm, the velocity
of the top roller is Vtop = 0.07366 m/s and for the bottom roller Vbottom = 0.07383 m/s. For
the point of attachment it was approximated, based also on the experimental data of
Vlachopoulos and Hrymak, that it is nine times larger than the point of detachment so it
was placed at xi = -37.31 mm (the half–thickness of the sheet at xi is Hi = 6.00 mm).

Figure 6.6. Top view of the computational domain for the spreading of the sheet.

58

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.7. Inclined–top view of the computational domain for the spreading of the sheet.

It is also assumed that the feeding of the material, is done from a single screw extruder
which is located far behind the calenders. Between the extruder and the calenders it is
assumed that a rectangular channel (i.e. a flat die) exists, which is attached directly to the
calenders. The width of the channel is Wen. and the gap is 2Hi. By this, it is assumed that
the material is confined before it enters to the region between the calenders. A rough
sketch of the above is shown at Figure 6.8. A similar calendering process is carried out by
Michelis [53] for the production of polymeric strips a few centimeters wide. In the
present case the simulation is performed for a Newtonian fluid with a relatively low
viscosity, that is μ=1000 Pa·s and with ρ=1000 kg/m3.

One of the most important steps for the simulation, after the construction of the
computational flow domain, is to decide the nature of the boundary conditions from a
physical point of view. The notation of each boundary can be found in Figure 6.9. In
reality, the side boundaries of the domain, where the spreading of the sheet occurs, are
free surfaces. A rational approximation that can be done is to accept that these boundaries
are in contact with the atmosphere, thus, the pressure can be set equal to zero. A similar
boundary condition for side boundaries for the pressure is followed by Luther et al [21]
and by Levine et al. [47].

59

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.8. Schematic representation of the process employed. The molten material is directly fed from the
extruder to the calenders through the rectangular channel. R is the radii of the calenders and U the velocity.

Moreover, it is also assumed that along the side boundaries, the fluid can slip. For that
reason the sides are approximated as surfaces with perfect slip. As stated by Weller [28,
29] the slip boundary condition for the velocity implies that the normal component has a
fixed value and equal to zero (no cross flow takes place) and that the tangential
components have zero gradients. Mathematically this is formulated as

U t
Un  0 and 0 (6.1)
n

Concerning the entrance, as in the two dimensional Newtonian case, the boundary
condition for the velocity is the fixed value of the flow rate. For the pressure boundary
condition it was decided to set the pressure gradient equal to zero. Mathematically, this
can expressed as
p
0 (6.2)
n

60

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.9. Notation of the boundaries on the computation flow domain.

Below, the reader can find summarized the above boundary conditions along with the
ones for the rest of the boundaries of Figure 6.9.

 Top and bottom calenders


p
- Pressure: Normal gradient of pressure is zero, 0
n
- Velocity: tangential velocity Ut ,rol.  V , normal velocity U n,rol.  0

 Sides
- Pressure: Fixed value, p = 0
U t
- Velocity: U n  0 and  0 (Un is the normal component and Ut the
n
tangential component)

 Outlet
- Pressure: Fixed value, p = 0
U
- Velocity: Normal gradient of velocity is zero, 0
n
61

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
 Inlet
p
- Pressure: Normal gradient of pressure is zero, 0
n
- Velocity: Known flow rate Q   Ut ,rol. Wen. 2H i 

where 2Hi is the distance between the rollers at the inlet.

6.1.2 Results and Discussion

The pressure field of the simulation as obtained with OpenFOAM is shown at


Figure 6.10. From a very first glance it can be seen that, at both side boundaries and near
the exit, the pressure exhibits negative values. This may be a result of the imposed
boundary conditions, combined with the geometry of the computational flow domain.
These regions, although they exhibit extremely negative values of pressure, with respect
to the maximum positive pressure developed, they are contained only in a small volume
of the computational flow domain something which is obvious from Figure 6.10. These
regions appear to occur near the exit boundary and adjacent to the two side boundaries.
The regions appear to be symmetric with respect to the centerline.

In the real calendering process, the sheet does not spread linearly from the
entrance to the exit. Such a claim holds valid, as in Unkrüer’s work [13] it is reported that
in the narrow region of the nip (minimum distance between the calenders), drag flow
caused by the movement of the calenders is predominant as compared to the flow in the
cross – machine direction, thus, it can be concluded most of the sheet spreading occurs
before the material reaches the nip region.

The pressure profiles at different locations are shown at Figure 6.11 and at Figure
6.12 plotted in a three–dimensional form. The pressure profile in the machine direction
appears to have the usual bell–shape behavior with the maximum pressure to be
approximately 3.17 MPa. It can also be seen that near the entrance the pressure exhibits a
build–up (short pressure plateau) of approximately 0.76 MPa as shown at Figure 6.11.

62

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.10. Calculated pressure profile for the linear geometry variation (width of sheet).

This is the result of the applied boundary condition of the pressure at the specific
boundary. Such a behavior is rational since it was assumed that the material is fed from a
rectangular die fitting directly to the calenders. The presence of the calenders and the
rectangular die, at the entrance, creates a resistance at the flow thus creating a pressure
build–up and this is exhibited by the plateau of the pressure curve, along the centerline
and near the entrance. After the built–up the pressure increases to a maximum value due
to the rotational movement of the calenders, and then drops to zero at the exit, as dictated
by the boundary condition on that boundary. Moreover, at Figure 6.11 the maximum
pressure is located at x=-xf =-4.37 mm, thus, validating once again the strength of the
Lubrication Approximation (as discussed in Chapter 2).

Centerline
3,5

3,0

2,5

2,0
Pressure (MPa)

1,5

1,0

0,5

0,0

-0,5
-0,04 -0,03 -0,02 -0,01 0,00 0,01
x - position (m)

Figure 6.11. Pressure distribution along the centerline located at the middle section of the computational
flow domain.

63

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
At Figure 6.12 the pressure distribution in the cross machine direction at different
locations is shown. The locations at which the pressure distribution was calculated in the
cross–machine direction were at x=0, x=-xf, x=-2xf and x=-3xf. The shape of the curve is
quasi–parabolic with a maximum at the middle section of the calenders (Wen./2=112.50
mm). As we move at the two sides, the pressure starts to level–off, becomes zero, then
takes negative values and finally increases to obtain again zero value. Pressure
distribution in the cross–machine direction for the position of x=0 (minimum distance
between the rollers), appears to have a more flatten profile than at other locations. This
effect is due to the rotational movement of the calenders, implying that the pressure
gradient does not change much for a large portion of the calenders, thus the predominant
flow is the one in the machine direction induced by the rotational movement. This
validates the claim presented previously that most of the material spreading occurs before
it reaches the nip region. However, there are many values of negative pressure at the
sides close to the nip region till the exit. On the other hand, pressure profiles which are
near the entrance exhibit a more quasi – parabolic profile and the values of the pressure
are all positive or negative but very close to zero.

Centerline
x=0.00 mm
-xf=-4.37 mm
-2xf=-8.74 mm
3
-2xf=-13.11 mm

2
)
Pres sure (M Pa

-1
0,01
0,00
m)

-2 -0,01
n(

-0,1
-0,02
t io

0,0
z-d
ec

0,1 -0,03
ir e c
d ir

t io n 0,2
(m )
x-

0,3

Figure 6.12. Pressure distribution along the x and z direction at different distances for the case of linear
width variation.

64

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
The pressure distribution in the cross–machine direction, at the nip region (x=0),
exhibits a quasi–parabolic profile which is in qualitative agreement with the results of the
pressure distribution along the axis of the calenders measured by Unkrüer as shown at
Figure 6.13 (The interest lies in the curve that belongs to ho = 0.25 mm since that half–
gap is very close to the one of the present study (i.e. ho = 0.30 mm)). The shape of the
curves in Unkrüer’s study is affected by the distance between the calenders and the
presence of the melt bank whose shape also alters in a certain amount the curves.
However, in the study presented it was assumed that near the entrance there is no melt
bank, the presence of which occurs in the experiment [11-13], which might affect the
results. Even though it was adopted that no melt bank was developed near the entrance
the results showed a remarkable qualitative agreement with Unkrüer’s experimental
work.

Figure 6.13. Pressure distribution in the cross – machine direction for different minimum gaps (ho)
between the calenders as measured by Unkrüer [13]. hk is the height of the melt bank and Vw1, Vw2 are the
velocity of the calenders.

From the same results, and since the objective of the thesis is to obtain the
pressure profile in the cross–machine direction, the imposed boundary conditions seemed
to be efficient to describe numerically the calendering process and produce some very
interesting and novel findings for the first time, using the fully three–dimensional

65

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Navier–Stokes equations. However, on the effort to predict also a possible shape for the
sheet spreading from the entrance to the exit of the flow domain, we need to take into
account the presence of the negative pressure regions. If we set as a point of reference
those regions and how to eliminate or to minimize them, it can be concluded that,
although the boundary conditions describe sufficiently the process, the linear spreading
seems to be only good as a starting guess. For those reasons, the next step is to construct
a flow domain at which the sheet spreads linear to a finite width from the entrance to
some point behind the nip region. From that point till the exit it is assumed that the
spreading is terminated and the width of the sheet remains the same until the exit.

Closer examination of the results of Figures 6.10 and 6.12, reveals that in the
machine direction, the negative pressure regions exist approximately between -2xf and xf.
An educated guess to eliminate those regions would be to construct the flow domain in
way that the width of the sheet geometry varies linearly from xi to -2xf and from -2xf
until the exit there is no any geometrical variation of the sheet’s width. The new flow
domain is shown at Figure 6.14.

Figure 6.14. Computational flow domain with the spreading of the sheet terminated before the exit.

For the computational flow domain constructed above the calculate pressure field
as obtained by the simulation is shown at Figure 6.15. The boundary conditions were
assumed to be the same as in the previous simulation. It is clear from the figure, that the

66

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
negative pressure regions which appeared in the previous simulation are minimized if
compared to Figure 6.10. These regions appear to be near the two side boundaries and at
the exit. The approximated shape of the sheet caused a decrease in the negative pressures
(absolute value) which means that the assumption of linear width variation to a specific
position before the nip region seems reasonable to be assumed. Of course this is only an
approximation that we assume and this is why regions with negative pressure still exist.
At the side boundaries, where the width of the sheet remains constant, no negative
pressure regions were observed. The negative values lie near the exit boundary and the
two sides where the width varies linearly.

Figure 6.15. Calculated pressure field for the case where the spreading terminates near the nip region.

The pressure distribution along different directions and positions is shown at


Figure 6.16. It can be seen that the plotted curves exhibit almost the same behavior as in
the previous analysis if compared to Figure 6.12 (the pressure distribution curves are
plotted at exactly the same positions as in the previous simulation concerning the linear
spreading from the entrance to the exit of the computational flow domain). Although the
shape of the flow domain was approximated with a different manner than the previous
numerical experiment, the pressure curves preserve their quasi – parabolic behavior (the
maximum pressure seems to be unchanged with regard to the previous simulation),
implying that the assumed boundary conditions describe satisfactory the problem. Even

67

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
closer comparison of Figure 6.12 and Figure 6.16, it can be seen at the latter figure, that
the pressure distribution in the cross–machine direction exhibits a “more” quasi–
parabolic shape, which means that the pressure gradients exhibit a larger change as we
move from the central region to the sides. Near the exit boundary, the pressure
distribution in the cross–machine direction appears to have a more flatten profile as it was
reported in the previous three–dimensional numerical experiment. Moreover, near the
sides the pressure values for the positions of -xf =-4.37 mm and x=0 are all positive. For
the position of -2xf =-8.74 mm, some of the calculated pressure values hold negative
numbers, the absolute values of which are smaller if compared to the maximum pressure
developed by the simulation.

Centerline
x=0.00 mm
3,5 -xf=-4.37 mm
3,0 -2xf=-8.74 mm

2,5
P a)

2,0
P re ss ure (M

1,5

1,0

0,5

0,0
0,000
-0,5
-0,2 -0,015
)
(m

-0,1
n

-0,030
io

z- 0,0
d ir e
ct

c t io 0,1
ir e

n -0,045
-d

0,2
x

Figure 6.16. Pressure distribution along the x and z direction at different distances for the case of linear
width variation from the entrance to -2xf.

The streamlines as calculated by the simulation are showed at Figure 6.17. The
behavior of the streamlines shows the spiral transport of the material from the entrance to
the two sides and to the exit of the computational flow domain. The material conveys
from the entrance of the domain to the two side boundaries and simultaneously
progresses to the exit of the domain. This spiraling motion is in qualitative agreement
with the numerical results by Luther and Mewes [21] and with the experimental results

68

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
observed by Unkrüer [13]. The only difference is that in the present analysis there is not
any melt bank. In a sense, it can be assumed that the spiraling motion in the present
simulation belongs to a hypothetical melt bank. Moreover, the spiraling is caused due to
the imposed boundary conditions for pressure and velocity at the side boundaries which
force the material to flow, except the machine direction, also in the cross–machine
direction implying the existence of velocity component uz. The velocity field for uz is
presented at Figure 6.18 for different cross–sectional areas normal to the z–axis and at
Figure 6.19 for a cross–section normal to the y–axis. It can be seen that uz is greater in
two regions, which appear to be symmetrical, located between the entrance and the two
sides. The negative values are due to the minus z–axis. In the middle region, uz is
approximately zero which means that the sheet does not spread close to that region and
this is why at Figure 6.17 the spiraling effect of the streamlines is less intense.

Figure 6.17. Calculated streamlines for the case where the spreading terminates at -2xf. The legend refers to
the magnitude of the velocity.

Another result of the imposed boundary conditions is that some of the streamlines
in Figure 6.17, terminate at the side boundaries while some others start from the entrance,
they progress tangentially to the side boundaries and they finally terminate at the exit
boundary of the domain. This implies that from the two sides there is outflow. The
outflow seems to be greater from the wider area of the side boundaries. From the
narrower area of the same boundary, the outflow is much less. This seems to be rational

69

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
since in the nip region the flow in the x–direction is predominant and the flow in the z–
direction is less intense. Of course, the assumptions done previously are only an initial
reasonable approximation to describe the flow between the calenders. The negativity of
the pressure seems to be more prone to the shape of the computational domain rather than
the imposed boundary conditions. However, these approximations can be refined,
especially for the geometry of the side boundaries and a parabolic or a quasi–parabolic
shape could probably eliminate or minimize the negative pressure regions. The present
study gave only a sufficient insight on how to construct an appropriate computational
flow domain to describe the hydrodynamics of the flow as well as the prediction of
pressure distribution in the cross–machine direction.

Figure 6.18. Top view of the velocity field for the velocity component uz for a cross – sectional area which
is normal to the y–axis.

70

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.19. Inclined view of the velocity field for the velocity component uz for different cross – sectional
areas normal to the z–axis.

In the present simulation, as well as in the previous one, the shape of the
computational domain was determined based on a set of boundary conditions and a pre–
determined spreading which was improved after a single simulation. The spreading of the
sheet is the result of the imposed boundary conditions combined with the shape of the
computational flow domain. The observed velocity in the cross–machine direction, uz, is
the result of the pre–determined geometry. Thus, the fluid must spread so that to fill the
expanding geometry.

6.2Numerical Results with Carreau Viscosity Model

6.2.1 Geometry, Spreading of the Sheet and Boundary Conditions

In the previous simulations, the Newtonian viscosity model gave very satisfactory
results, most of which are in qualitative agreement with past theoretical and experimental
studies [13,21,44]. In the present section the Carreau viscosity model is employed for the
simulations. The geometry for the computational domain is obtained from the study by
Vlachopoulos and Hrymak [14] and the spreading is considered to be the one presented in
the final three dimensional simulation above, for the Newtonian fluid where the

71

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
spreading terminates at -2xf. Thus, the computational domain is the one shown at Figure
6.14. The parameters for the Carreau viscosity model are the ones presented in the two–
dimensional study in Chapter 5. For the boundary conditions, it seemed rational to
impose the ones used in the three–dimensional Newtonian analysis of the previous
section. The way the material is fed in the inlet is the one described in the previous
analysis, illustrated at Figure 6.8. The aim of the present simulation is not to validate the
work done by Vlachopoulos et al. [14]. From their work we only obtained geometrical
and rheological parameters. Our formulation is in the scope of presenting qualitative
results for the pressure distribution in the cross – machine direction, spiraling movement
of the material near the entrance (with no presence of melt bank) and approximate
spreading of the sheet. Results by Unkrüer [13] and Luther et al. [21] provide the
necessary information for such a qualitative comparison.

6.2.2 Results and Discussion

The pressure distribution for the computational domain is presented at Figure


6.20. The computed field appears to have similar pressure behavior with the previous
three–dimensional simulations for the Newtonian fluid with the pressure now being
higher. It can be seen, in the same figure, the presence of negative regions as in the
previous simulations. Such a result was anticipated due to the approximations of the
boundary conditions and the shape of the computational flow domain, as it was discussed
in the previous section. For that reason, we can conclude that the viscosity model used
each time does not affect the presence of negative regions in the flow domain.

The pressure distribution in the machine and cross – machine direction at different
positions is presented in Figure 6.21. The maximum pressure is reported to be
approximately pmax=9.1 MPa. The pressure distribution along the centerline has the well
known bell–shape shape. Near the entrance there is a relatively large pressure built–up
which is caused by the obstruction of the flow, as it was assumed that exactly before the
calenders there is attached a rectangular die from where the polymer melt is fed. In the
nip region (x=0) the pressure profile in the cross – machine direction is more flatten if
compared to the rest of the curves of Figure 6.21, implying that the dominant flow is

72

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Figure 6.20. Calculated pressure field with the Carreau viscosity model.

due to the rotational movement of the calenders, thus the spreading of the sheet is small
as the pressure gradient does not change much in the z–direction. At x=-xf and for x<-xf
the pressure profiles exhibit a “more” quasi – parabolic profile, and the rapid change of
the pressure gradients imply the rapid spreading of the sheet. However, as we move near

10 Centerline
x = 0.00 mm
x = -4.37 mm
8 x = -8.36 mm
x = -13.11 mm
6
Pa )
Pres sure (M

0,01
0 0,00
m)

-0,01
n(

-0,15
t io

-0,10 -0,02
ec

-0,05
0,00 -0,03
d ir

0,05
z-d 0,10 -0,04
x-

ir e c 0,15
t io n
(m )

Figure 6.21. Pressure distribution in the machine and cross – machine direction at different positions with
the Carreau viscosity model.

73

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
the entrance, there are some negative pressure values near the two side boundaries. For
the present simulation it was accepted that the obtained absolute value of the negative
pressure is the maximum that can be numerically obtained with the present set of
boundary conditions and the computational domain geometry. The total elimination of
the negative regions would require a more thorough examination of the shape of the
computational domain, especially for the two side boundaries, and also decide where the
spreading of the sheet terminates. Of course, the lack of melt bank near the entrance of
the computational domain affects somehow the results but since the present analysis is on
the scope of a qualitative study, if someone compares Unkrüer’s results shown at Figure
6.13 a lot of qualitative similarities can be found with the present work. So far, there are
not any reports in the open literature for quasi – parabolic pressure profile in the cross–
machine direction using the fully 3D Navier – Stokes equations along with a Generalized
Newtonian Fluid model. Thus, we are the first to obtain qualitative results by the use of
the above approximations. A more elegant approach, could employ the determination of
the shape of the side boundaries on the basis that the width of the sheet does not vary
linearly, but in a parabolic or quasi – parabolic manner. The exact shape for the spreading
can be determined through several iterations concerning the shape of the side boundaries
in a way that the streamlines near the boundary do not cross it but lie tangentially to it.

Figure 6.22. Calculated streamlines for a polymer melt obeying the Carreau viscosity model.

74

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Furthermore, at Figure 6.22, the calculated streamlines of the flow domain is
presented. In the central region the spiraling effect is less intense and it becomes greater
as the material progresses to the two side boundaries. Some of the streamlines end
directly to the side boundaries and some others end at the exit boundary. For those
streamlines that terminate at the side boundaries, the material can outflow from the
greater area of the sides, whereas from the smaller area the outflow, if any, is small due to
the fact that the rotational movement of the calenders create a dominant drag flow in the
nip region that forces the material to flow in the machine direction rather than spread.

For those reasons, the side boundaries are not treated, strictly speaking, as free
surfaces. The partly elastic nature polymer melts hold, prevents them to spread infinitely
on such a process. Thus, the present boundary conditions seemed rational to be assumed.
In a matter of speaking it is approximated that the liquid–elastic material is treated, with a
purely viscous model like most fluids. To account for the finite spreading we assumed a
set of boundary conditions that seemed to be rational from a physical point of view, and
imposed them on the spirit of preventing the material to spread infinitely. Of course,
outflow from the side boundaries is not a desirable parameter. However, the results
showed qualitative agreement with theoretical and experimental past studies but there is
always room for improvement.

As a final result of the present numerical experiment, we show the velocity


profiles of Umag., ux and uz at different locations of the computational domain, along the
y–direction. The different locations are represented schematically at Figure 6.23.

Figure 6.23. Schematic representation of the locations for the velocity profiles along the y–axis.

75

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Umag.
0,130 ux
uz

0,128

y - position (m)
0,126

0,124

0,122

0,120
-0,01 0,00 0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08
Velocity (m/s)

Figure 6.24. Velocity profiles for location (A).

At Figure 6.24, ux is higher near the calenders and decreases in the mid–plane (the
plane normal to y–axis) taking negative values. The negative values are due to the
recirculation patterns predicted by the numerical experiment and imply backflow in the
machine direction. This is also why Umag. profile has the shape shown at the figure which
means that the material’s velocity increases to a certain value near the centerline but it
moves backwards. Simultaneously, at that location (A), uz holds zero (or very small)
velocity values which means that the material does not flow in the cross machine
direction, thus does not spread. If we move from that location a few centimeters, either to
the one side or the other, uz holds higher values which means the sheet is assumed to
spread thus resulting in the spiral motion discussed above and showed at Figure 6.22
which appears to be weak close to location (A) and intense as the material progresses to
the sides.

In Figure 6.25, at position (B), the velocity component uz, appears to have a
quasi–parabolic pressure profile similar to the well known parabolic velocity profile for a
pressure driven flow between two flat plates. The velocity magnitude Umag. profile in the
same figure, seems to have higher values than in location (A) and this is due to the
contribution of the velocity in the z–direction as the fluid is conveyed downstream and in

76

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
the machine direction. This velocity, as discussed extensively in previous sections, is the
result of the applied boundary conditions at the side boundaries of the flow domain and

Umag.
0,130 ux
uz

0,128
y - position (m)

0,126

0,124

0,122

0,120
-0,01 0,00 0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08
Velocity (m/s)

Figure 6.25. Velocity profiles for location (B).

Umag.
0,1275 ux

0,1270
uz

0,1265

0,1260
y - position (m)

0,1255

0,1250

0,1245

0,1240

0,1235

0,00 0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08


Velocity (m/s)

Figure 6.26. Velocity profiles for location (C).

the shape of the geometry. Similar behaviors seems to have the velocity profiles at
location (C) as shown in Figure 6.26. The velocity component uz, appears to have a
quasi–parabolic profile which means that the material is conveyed in the z–direction.
From the ux velocity profile, it can be seen that there are no negative values which means
that there is not any back flow, thus the circulatory patterns do not occur at that location.

77

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Umag.
ux
0,1257
uz
0,1256

0,1255

0,1254

y - position (m)
0,1253

0,1252

0,1251

0,1250

0,1249

0,00 0,01 0,02 0,03 0,04 0,05 0,06 0,07


Vel Mag

Figure 6.27. Velocity profiles for location (D).

The velocity magnitude Umag. profile shown at the same figure, reveals that the velocity is
maximum near the calenders and drops to an approximately constant value near the mid–
plane x–z. Near that plane, the velocity is almost the same to the maximum velocity in
the z–direction, thus the dominant flow occurs in that direction.

The final plotted velocity profiles are near the exit, at location (D), and are
presented at Figure 6.27 with the velocity profile downstream being almost constant
(Umag.=ux=0.0783 m/s) and the velocity in the cross – machine direction very low or zero
values. This implies, that the sheet at the exit has very close velocity with the velocity of
the calenders.

78

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Closure and Future Work

The proposed objectives of the present thesis were to predict the pressure
distribution in the cross–machine direction as well as a sufficient geometry with varying
width of the calendered sheet by the help of OpenFOAM. The study was initiated by
validating a previous work on calendering using the Newtonian and the Generalized
Newtonian Fluid models. For the Newtonian Fluids the numerical experiments exhibited
excellent agreement with the Lubrication Approximation Theory. However, the
peculiarity of the power–law model, combined with the lack of measured rheological
data, gave erroneous results and for that reason the Carreau viscosity model was utilized
giving results with good overall agreement. The study for the three–dimensional flow
between the calenders for Generalized Newtonian Fluids, with the assumption that no
melt bank occurs, gave results which were in qualitative agreement with past numerical
and experimental studies. The results gave some very interesting insight into the quasi–
parabolic pressure distribution in the cross–machine direction and the shape of the
calendered sheet as it spreads. Thus, the present work is so far the only one that can
predict the above, employing the Navier – Stokes equations by means of a Generalized
Newtonian Fluid model and making several reasonable approximations throughout the
study. Of course, there is always room for improvement.

The inability of the Navier–Stokes equations to predict the partly elastic behavior
of the material is a limiting parameter for such a study. However, by making special
assumptions the material can be treated as a purely viscous fluid. In future work, special
care needs to be taken for the side boundaries. The linear width variation proposed here,
is only the one side of the coin. The other side would require an effort to express those

79

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
sides with a parabolic or quasi–parabolic shape. This can be performed by a certain
number of simulations in the scope of the calculated streamlines do not cross the side
boundaries, thus no outflow from these sides should occur. Moreover, another
assumption that needs to be taken into account is the presence of melt bank and on what
extend its shape affects the results. Finally, another important aspect that someone should
embody in the study is the non–isothermal study of the process. This can be done by
inserting and solving simultaneously with the Navier–Stokes, the energy equation and
also introducing temperature dependence of the viscosity models.

80

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
Bibliography

1. Middleman S., Fundamentals of Polymer Processing, McGraw–Hill, New York,


USA, 1977
2. Tadmor Z., Gogos C., Principles of Polymer Processing, John Wiley & Sons,
New Jersey, USA, 2006.
3. Gujati D.P., Leonov A.I., (Eds), Modeling and Simulation in Polymers, Wiley–
VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2010.
4. Cogswell F.N., Melt Rheology, Woodhead Publishing, Cambridge, 1996
5. Gaskell R.E., The Calendering of Plastic Materials, J. App. Mech.,17, 334–337,
1950
6. McKelvey J.M., Polymer Processing, Wiley, New York, USA, 1962.
7. Brazinsky I., Cosway H.F., Valle Jr. C.F., Clark R., Jones R., Story V., A
Theoretical Study of Liquid–Film Spread Heights in the Calendering of
Newtonian and Power–Law Fluids, J. App. Polym. Sci., 14, 2771–2784, 1970
8. Alston W.W., Astill K.N., An Analysis for the Calendering of Non–Newtonian
Fluids, J. App. Polym. Sci., 17, 3157–3174, 1973
9. Kiparissides C., Vlachopoulos J., Finite Element Analysis of Calendering,
Polymer Engineering and Science, 16 (10), 712–719, 1976
10. Kiparissides C., Vlachopoulos J., A Study of Viscous Dissipation in the
Calendering of Power–Law Fluids, Pol. Eng. Sc., 18 (3), 210–214, 1978.
11. Agassant J.F., Espy M., Theoretical and Experimental Study of the Molten
Polymer Flow in the Calender Bank, Pol. Eng. Sci. 25 (2), 118–121, 1985.
12. Agassant J.F., Le Calandrage des Matières Thermoplastiques, Doctoral Thesis,
Paris 6, France, 1980.
13. Unkrüer, W., Beitrag zur Ermittlung des Drucksverlaufes und der Fliessvorgange
im Walzspalt bei der Kalanderverarbeitung von PVC Hart zu Folien, PhD Thesis,
IKV, TU Aachen, 1970.
14. Vlachopoulos J., Hrymak A.N., Calendering Poly(Vinyl Chloride): Theory and
Experiments, Pol. Eng. Sci., 20 (11), 725–731, 1980.

81

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
15. Mitsoulis E., Vlachopoulos J., and Mirza, F.A., Calendering Analysis Without the
Lubrication Approximation, Polym. Eng. Sci., 25, 6–18, 1985.
16. Vlachopoulos, J. and Mitsoulis, E., Fluid Flow and Heat Transfer inCalendering:
A Review, Transport Phenomena in Polymeric Systems–2, vol. VI (eds. M.R.
Kamal, R.A. Mashelkar, and A.S.Mujumdar), Advances in Transport Processes,
Wiley Eastern, New Delhi, 79–104, 1988.
17. Sofou S., Mitsoulis E., Calendering of Pseudoplastic and Viscoplastic Sheets of
Finite Thickness, Journal of Plastic Film & Sheeting, 20 (3), 185 – 122, 2004
18. Mitsoulis E., Sofou S., Calendering Pseudoplastic and Viscoplastic Fluids with
Slip at the Roll Surface, J. App. Mech., 73 (2), 291–299, 2006
19. Kovacs J., Mandjarov R., Blisnjuk T., Prehn K., Sussiek M, Muller J, Schulte K.,
Bauhofer W, On the Influence of Nanotube Properties, Processing Conditions
and Shear Forces on the Electrical Conductivity of Carbon Nanotube Epoxy
Composites, 20,1–6, 2009.
20. Syehan T., Gojny F., Tanoglu M., Schulte K., Critical Aspects Related to
Processing of Carbon Nanotube/Unsaturated Thermoset Polyester
Nanocomposites, Eur. Polym. J., 43(2), 374–379, 2007.
21. Luther S., Mewes D., Three–Dimensional Polymer Flow in the Calender Bank,
Pol. Eng. Sc. 44 (9),1642–1647, 2004.
22. Patankar, S.V., Spalding D.B., A Calculation Procedure for Heat, Mass, and
Momentum Transfer in Three–Dimensional Parabolic Flows, Int. J. Heat Mass
Transfer, 15, 1787–1806, 1972.
23. Weller H., OpenFOAM User Guide, OpenCFD Limited, www.opencfd.co.uk,
2004.
24. Jasak H., Error analysis and estimation for the finite volume method with
applications to fluid flows, PhD Thesis, Imperial College London, 1996.
25. Versteeg H.K., Malalasekera W., An Introduction to Computational Fluid
Dynamics The Finite Volume Method, Longman Scientific & Technical, London,
UK, 1995
26. Hestenes M.R., Stiefel E., Methods of Conjugate Gradients for Solving Linear
Systems, J. Res. Natl. Bur. Stand., 49, 409–436, 1952

82

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
27. Breuer M., Direkte Numerische Simulation und Large–Eddy Simulation
turbulenter Strömungen auf Hochleistrungsrechnern, Habilitationsschrift,
Universität Erlangern–Nuerberg, Berichte aus der Stroemungstechnik, Shaker
Verlag, Aachen, 2002.
28. Biswas G., Breuer M., Durst F., Backward–Facing Step Flows for Various
Expansion Ratios at Low and Moderate Reynolds Numbers, J. Fluids Eng., 126,
362–374, 2004.
29. CALENDERCAD software, Polydynamics, Inc., Dundas, Ontario, Canada
30. Weller H., OpenFOAM Programmer’s Guide, First Edition, OpenCFD Limited,
www.opencfd.co.uk, 2004.
31. Zheng R., Tanner R.I., A Numerical Analysis of Calendering, J. Non–Newtonian
Fluid Mech., 28, 149–170, 1988.
32. Mitsoulis E., Finite Element Analysis of Two–Dimensional Polymer Melt Flows,
PhD Thesis, McMaster University, 1984.
33. Mitsoulis E., Numerical Simulation of Calendering Viscoplastic Fluids, J. Non –
Newtonian Fluid Mech., 154, 77–88, 2008.
34. Silliman W.J., Scriven L.E., Separating Flow Near a Static Contact Line: Slip at
a Wall and Shape of a Free Surface, J. Comp. Phys., 34, 287–313, 1980.
35. Jay P., Piau J.M., Kissi N. E., Cizeron J., The Reduction of Viscous Extrusion
Stresses and Extrudate Swell Computation Using Slippery Exit Surfaces, J. Non–
Newtonian Fluid Mech. 79, 599–617, 1998.
36. Mitsoulis E., Some Issues Arising in Finding the Detachment Point in
Calendering of Plastic Sheets, J. Plastic Film and Sheeting, 26, 141–165, 2010.
37. Peck M.C., Rough S.L., Wilson D.I., Roller Extrusion of a Ceramic Paste, Ind.
Eng. Chem., 44, 4099–4111, 2005.
38. Kim N, Kobayashi S., Three dimensional simulation of gap controlled plate
rolling by the finite element method. Int. J. Mach. Tools Manuf., 30, 269–81,
1990.
39. Hwang S.M., Joun M.S., Kang Y.H., Finite element analysis of temperatures,
metal flow and roll pressure in hot strip rolling, Trans. ASME, J Ind., 115, 290–
98, 1993.

83

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7
40. Dvorkin E.N., Goldschmit M.B., Cavaliere M.A., Amenta P.M., Marini O.,
Stroppiana W., 2D finite element parametric studies of the flat rolling process., J.
Mater. Process Tech., 99–107, 1997.
41. Lee J.D., A large strain elastic–plastic finite element analysis of rolling process,
Comp. Met. App. Mech. Eng., 161, 315–47, 1998.
42. Komori K., Rigid-plastic finite element method for analysis of three dimensional
rolling that requires small memory capacities, Int. J. Mech. Sci., 40, 479–91,
1998.
43. Jiang Z.Y., Tieu A., A simulation of three-dimensional metal rolling processes by
rigid-plastic finite element method, J. Mater. Process Technol., 112, 144–51,
2001.
44. Duan X., Sheppart T., Three dimensional thermal mechanical coupled simulation
during hot rolling of aluminium alloy 3033, Int. J. Mech. Sci., 44, 2155–72, 2002.
45. Jiang Z.Y., Tieu A.K., Zhang X.M., Lu C., Sun W.H., Finite element simulation
of cold rolling of thin strip., J. Mater. Process Technol., 140, 542–7, 2003.
46. Matsumiya T., Flemings M.C., Modeling of Continuous Strip Production by
Rheocasting, Metallurgical Transactions B, 12B, 17–31, 1981.
47. Sezek S., Aksakal B., and Can Y., Analysis of Cold and Hot Plate Rolling Using
Dual Stream Functions, Materials and Design, 29, 584–596, 2008.
48. Levine L., Corvalan C.M., Campanella O.H., Okos M.R., A Model Describing the
Two–Dimensional Calendering of Finite Width Sheets, Chem. Eng. Sci., 57, 643–
650, 2002.
49. Teletzke G.F., Davis H.T., Scriven L.E., How Liquids Spread on Solids, Chem.
Eng. Comm., 55, 41–81, 1987.
50. Huh C., Scriven L.E., Hydrodynamic Model of Steady Movement of a
Solid/Liquid/Fluid Contact Line, 35, 1, 85–101, 1971.
51. Li D., Gu Y., A model for liquid drop spreading on a solid surface, Colloids
Surfaces A: Physicochem. Eng. Aspects, 142, 243–256, 1998.
52. Ferry J.D., Viscoelastic Properties of Polymers, 3rd Edition, John Wiley & Sons,
U.S.A.,
53. Michelis P., IMMG S.A., Personal communication.

84

Institutional Repository - Library & Information Centre - University of Thessaly


08/12/2017 18:41:56 EET - 137.108.70.7

You might also like