Three-Dimensional Flow Analysis in The Calendering Process
Three-Dimensional Flow Analysis in The Calendering Process
Three-Dimensional Flow Analysis in The Calendering Process
School of Engineering
Department of Mechanical Engineering
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
ii
iii
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
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
vii
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].
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].
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 kT 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
U 0 (2.4)
0 p (2.5)
C p U T k2T : 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)
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
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
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
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
where
1 u u 1 u u
Dxx x x , Dxy x y (2.20)
2 x x 2 y x
2 D (2.21)
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
mT
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
1
ΙΙ D (2.24)
2
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
For α=2 the model is simply called Carreau. The Cross model is the following
(2.26)
1
1n
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
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 hx and 0 on y 0 (2.31)
y
11
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
hx 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
hx 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)
12
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
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 32 1 52
x '
1 32
tan 1
x ' tan
1
1 32
(2.39)
32 H o 1 x' 22
1 2
13
Figure 2.2. Schematic representation of the pressure distribution between the calenders as sketched by
Middleman [1].
3C R
pmax ' (2.40)
2 2H o
where
1 32
C
1 2
1 32 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 32 2 1
p' x
9R
1 3 tan (2.42)
2
32 H o 1 2
14
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
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.
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 hx 1 n (2.48)
q m dx
16
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
n
pH
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.
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
17
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].
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
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].
19
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
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
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
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.
In all the studies presented above, the calendering process was studied using two–
dimensional numerical schemes (mostly Finite Elements Method) or one–dimensional
23
24
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].
25
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.
26
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)
27
C
C C
C C (4.5)
C C
C C
C
1
2 (4.6)
28
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.
The aim is to minimize the value of the residual. If that value is small enough, then it can
be approximated that
29
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].
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.
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
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.
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.
31
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].
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
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.
33
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
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
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.
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
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
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
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).
39
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
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
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
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
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
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
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
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
Figure 5.14. Velocity components profiles and velocity magnitude at x=6.72 mm (exit boundary).
46
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
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
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
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
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].
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
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
52
3D Effects in Calendering
Introduction
53
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
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
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].
56
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
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
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
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
Below, the reader can find summarized the above boundary conditions along with the
ones for the rest of the boundaries of Figure 6.9.
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
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
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
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
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
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
Figure 6.15. Calculated pressure field for the case where the spreading terminates near the nip region.
67
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
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
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
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.
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
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
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
Figure 6.22. Calculated streamlines for a polymer melt obeying the Carreau viscosity model.
74
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.
Figure 6.23. Schematic representation of the locations for the velocity profiles along the y–axis.
75
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)
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
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)
Umag.
0,1275 ux
0,1270
uz
0,1265
0,1260
y - position (m)
0,1255
0,1250
0,1245
0,1240
0,1235
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
0,1255
0,1254
y - position (m)
0,1253
0,1252
0,1251
0,1250
0,1249
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
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
80
81
82
83
84