Viscoelastic Simulation of Bi-layer Coextrusion in a Square Die: an Analysis of Viscous Encapsulation
Mahesh Gupta
Michigan Technological University Plastic Flow, LLC
Houghton, MI 49931 Hancock, MI 49930Abstract
The Giesekus model is used for viscoelastic
simulation of a bi-layer flow in a square die. In contrast to
the experimental data reported in the literature, in the
present work even with viscoelastic effects included in the
simulation, encapsulation of a high viscosity polymer by a
lower viscosity polymer could not be captured. Since the
viscous encapsulation could not be captured with a purely
viscous formulation either, it is concluded that the
difference in the wettability and surface tension of the two
polymers is probably the major factor resulting in the
encapsulation.
Introduction
Coextrusion, which involves simultaneous extrusion
of several different polymers through a die, combines the
functionalities and benefits of several polymers into a
single multi-layered product [1]. Despite this inherent
advantage of coextrusion, design of coextrusion dies is
difficult because depending upon the rheology of
polymers used for coextrusion, the polymers in various
layers may get redistributed as they flow through the die
such that the distribution of various polymers at the inlet
and at the exit of the die may be quite different. In
particular, it has been observed in coextrusion
experiments that as two polymers with different
viscosities flow together, in the region near the die walls
the lower viscosity polymer tends to encapsulate the
higher viscosity polymer [2]. Because of the undesirable
variation in layer thickness caused by the encapsulation,
plastic films and sheets are often trimmed at the ends to
discard these non-uniform portions.
Unfortunately, the root cause behind the encapsulation
phenomena, whether it is caused by the difference in
viscosity, or in viscoelasticity, or in some other property
of the two polymers, is still not completely understood.
Simulations of coextrusion using a purely viscous
formulation available in the literature [3], including our
earlier work [4], have not been successful in capturing the
encapsulation. This may lead one to believe that the
encapsulation must be caused by viscoelastic effects,
which was supported by some papers [5] on viscoelastic
simulation of a bi-layer flow in the literature. In contrast,
the simulation results presented in this paper indicate that
the encapsulation is not caused by viscoelastic effects.
Therefore, at this point we believe that the unbalanced
force at the contact line where the two polymers together meet the die wall is the main driving force which results
in the observed encapsulation. This unbalanced force at
the contact line originates from the difference in the
wettability and surface tension of the two polymers.
Viscoelastic Formulation for Polymeric Flow
To capture the viscoelastic behavior of polymers, the
Giesekus model (Eqn. 1) [6] was employed in the current
work.
Besides the constitutive equation, simulation of
viscoelastic flow requires solving the mass and
momentum conservation equations. Assuming a steady,
inertia-less, isothermal, incompressible flow with no body
force, the conservation equations for momentum and mass
are simplified to
where the total stress, σ is given by the following
equation
with p being the pressure, n
s
the solvent viscosity, δ is the
identity tensor, and is the number of modes in the
viscoelastic constitutive equation.
Numerical Scheme for Flow Simulation
To obtain the finite element formulation of the flow
problem, in the present work, the standard Galerkin
method [7] was used for mass and momentum
conservation equations (Eqns. 3 – 5). However, for
stability of the numerical scheme, a streamline-upwind- Petrov-Galerkin (SUPG) method [8] was employed for the
Giesekus constitutive equation (Eqn. 1). For additional
stability of the numerical scheme the DEVSS method [9]
was employed for discretization of the momentum
equation (Eqns. 3 and 5). In DEVSS method an elliptic
operator 2η
a(e-e
a) is introduced in the momentum
equation, where e
a is a discrete approximation of the
strain-rate tensor e, and η
a = (sum) N n=i η
i
can be viewed as an artificial viscosity with η
i being the same as that in
Eqn. (1). For the exact solution to the flow problem this
additional elliptic operator is a null operator. However, in
the DEVSS method, e
a is explicitly determined by
solving the following equation.
In a finite element simulation of viscoelastic flow, strainrate
tensor is discontinuous across the element boundaries,
whereas a continuous interpolation is used for e
a,
resulting in a non-zero value of the additional elliptic
operator. In the present work, linear continuous
interpolation over tetrahedral finite elements was used for e
a. The linear interpolation over tetrahedral finite
elements was also used for pressure and the extra stress
tensor τ
i. However, for velocity, the linear continuous
interpolation over tetrahedral finite elements was
augmented by an extra node at the centroid of the
tetrahedral finite element. The extra node at the centroid
of each tetrahedral finite element is required to satisfy the
Babuska-Brezzi stability condition [10] for the
incompressible flow simulation.
Coextrusion Simulation Equations
For simulation of a multilayer flow of polymers
during coextrusion, the velocity and stresses are required
to be continuous across the interface between the adjacent
polymer layers [3]. That is,
u
(1) and τ
(1) are the velocities and stresses on one side of
the interface, with u
(2) and τ
(2) being the velocities and stresses on the other side of the interface.
Besides the continuity of velocity and stress, coextrusion simulation
requires enforcement of the no-cross-flow condition at the
interface. That is, the velocity component normal to the
interface must be zero at every point on the interface.
where u is the velocity, and n is the unit vector
perpendicular to the interface.
Mesh Partitioning Technique
In the three-dimensional simulations of coextrusion
reported in the literature, finite element mesh is modified
after each flow simulation iteration, such that the interelement
boundaries coincide with the interface between
adjacent layers of different polymers [3]. Such an approach using an interface-matched finite element mesh
can only be employed for simulating a two-dimensional
system or a simple three-dimensional system such as a
rectangular die. For real-life coextrusion systems, with
complex three-dimensional die channel geometry,
repeated generation and modification of interface-matched
finite element meshes is impractical.
In the present work, a three-dimensional mesh of
tetrahedral finite elements was generated over the
complete flow channel in the die. This finite element
mesh is not modified or regenerated at any stage during
coextrusion simulation. Thereby, allowing simulation of
even highly complex coextrusion systems.
In the mesh partitioning technique which is employed
in this work, the interface between adjacent layers of
different polymers is represented by a surface mesh of
linear triangular finite elements. However, the surface
mesh of triangular elements on the interface and the threedimensional
mesh of tetrahedral elements in the
coextrusion die are completely independent of each other.
This decoupling between the two finite-element meshes is
possible because in the mesh partitioning technique for
coextrusion simulation, the interface between adjacent
polymer layers is not required to match with the interelement
boundaries in the three-dimensional mesh of
tetrahedral finite elements. Instead, in the software used in
this work, the interface is allowed to pass through the
interior of the tetrahedral finite elements in the threedimensional
mesh.
In the mesh partitioning technique for coextrusion
simulation the tetrahedral elements which are intersected
by the mesh of triangular elements on the interface are
partitioned into two tetrahedral, pyramidal, or prismatic
finite elements. Further details of the mesh partitioning
technique are available in our earlier publications [4, 11].
Materials
For the bi-layer flow presented in this paper, Giesekus
model parameters given below were used for the two
different grades of polystyrenes [5].
The experimental data for the viscosities of the two
polystyrenes [3], along with the fit to the experimental
data using the Giesekus model parameters given above
[5], is shown in Fig. 1. It is evident from Fig. 1 that the
viscosity of Styron 472 is much higher than that of Styron
678E.
Bi-layer Coextrusion in a Square Channel
One of the main motivations for including viscoelastic
effects in coextrusion simulation in this work was to be able to capture encapsulation of one polymer by another
polymer which is often observed in polymer coextrusion.
For instance, Karagiannis et al. [3] obtained the
experimental data on encapsulation of Styron 472 by
Styron 678E in a bi-layer flow in a square channel.
Karagiannis at el. [3] observed that starting with a straight
interface shape at the contact line, where the two
polymers meet for the first time, Styron 678E
progressively encapsulated Styron 472 as the two
polymers flowed side-by-side in the square channel (Fig.
2). However, no such encapsulation was obtained in their
simulation of the bi-layer flow using a purely viscous
generalized Newtonian formulation (Fig. 3). Similar
simulation of bi-layer flow in a square channel using a
purely viscous formulation performed in our earlier work
[4] could not capture the encapsulation observed in
experiments.
It was reported by Sunwoo et al. [5] that they could
capture the encapsulation of Styron 472 by Styron 678E
when they included viscoelastic effects in the simulation,
which was the main motivation for including the
viscoelastic effects in coextrusion simulation in this work.
Sunwoo et al. [5] reported that the predicted encapsulation
in their simulation of the bi-layer flow was the largest and
the closest to the experiments when a large value of the
parameter α
i = 0.4 in Eqn. 1 was used in the simulation.
Therefore, in the present work, the same Giesekus model
parameters with α = 0.4, as those employed by Sunwoo
et al. for Styron 472 and Styron 678E (given above) were
also used in the present work to simulate the bi-layer flow
in the square channel. The simulation results thus obtained
are presented next in this section.
To characterize the flow, a non-dimensional shear rate
(Deborah number, De) is defined as
where U is the average velocity in the square portion of
the flow channel, b is the half of the square cross-section
side length, and λ
0 is characteristic relaxation time of the
polymer at low shear rate.
To define the Deborah number for the bi-layer flow in the
square channel, the average value of the characteristic
relaxation time, λ
0, for the two fluids employed was used.
For the average velocity in the square channel required for
a specified Deborah number, the appropriate value of the
uniform velocity was specified at the entrances of the two
polymers, whereas all components of the extra stress
tensor were specified to be zero at the entrances. Zero
velocity was enforced at the die walls.
For De = 1, with Styron 472 in the lower layer and
Styron 678E in the upper layer, the axial velocity
distribution for the bi-layer flow in the square channel is
shown in Fig. 4. Beyond De = 1 simulation of the bilayer
flow did not converge. Starting with the same
uniform velocity at the entrances of the two polymers,
after the two polymers meet the velocity slowly increases as the two polymers go through the converging portion of
the channel, and the highest velocity is obtained once the
two polymers reach the square portion of the channel. The
velocity then remains relatively unchanged as the two
polymers flow through the square channel
The transverse velocity distribution due to secondary
flow in the bi-layer coextrusion in the square channel for
De = 1 is shown in Fig. 5. For the bi-layer flow, in Fig. 5
there are eight recirculating vortices. However, in contrast
to the eight vortices in flow of a single polymer, which
have the same size, the vortices for the bi-layer flow in
Fig. 5 have different sizes. For instance, in Fig. 5 the two
vortices at the bottom are significantly smaller than the
two vortices just above the two bottom vortices. Also, in
contrast to the vortices for flow of a single polymer, for
the secondary flow in Fig. 5 there is no stagnation region
in the middle of the square cross-section. Instead, near the
center of the square cross-section the fluid from the
bottom vortices is entering into the flow in the two
vortices on top.
The shape of the interface between the two polymer
layers for De = 1 is shown in Fig. 6 (a). The
corresponding final shape of interface at the exit of the
square channel is shown in Fig. 6 (c). Starting with a
straight line interface shape at the contact line, where the
two polymers meet for the first time, the interface in Fig.
6 (a) starts to wave upwards near the middle and near the
two ends with troughs in between, resulting in the W-shaped
interface at the die exit in Fig. 6 (c). In contrast to
the wavy interface shape in Fig. 6 (a) and 6 (c), for De = 0.1 in Figs. 6 (b) and 6 (d), the predicted interface
shape just has a slight curvature but does not have a wavy
shape. For De = 0.1 the elastic effects in the flow, and
hence, the secondary flow vortices are negligible,
resulting in the simple shape of the interface which is
similar to the interface shape obtained in a purely viscous
simulation by Karagiannis et al. [3].
The magnitude of the second normal stress difference, |τ
22-τ
33|,for the bi-layer flow for De = 1 is shown in
Fig. 7. The magnitude of the second normal stress
difference, which results in the secondary flow vortices
observed in Fig. 4, is the maximum near the middle of
each of the four sides of square cross-section, and is zero
near the center of the cross-section, and is also zero along
the two diagonals of the square cross-section.
The pressure variation for the bi-layer flow in square
channel is shown in Fig. 8. As expected, the pressure is
zero at the exit and increases towards the two entrances.
Also, for the same velocity at the entrances of the two
polymers, the pressure gradient is higher in the narrower
feed channel of the upper layer.
Discussion
It should be noted that in contrast to the results
reported by Sunwoo at el. [5], no encapsulation of Styron
472 by Styron 678E was obtained in Figs. 6 (a) and (c).
For encapsulation, in Figs. 6 (a) and (c) the interface between the two layers should have bent downwards near
the two walls. Instead, in Figs. 6 (a) and (c), near the two
side walls the interface is bent upward, which is caused by
the secondary flow vortices shown in Fig. 4. Sunwoo at el.
[5] did not show the structure of secondary flows in their
simulation. Therefore, it is difficult to interpret how
Sunwoo et al. [5] could have obtained the encapsulation in
the bi-layer flow using the Giesekus model parameters
given in Section 3. Since encapsulation was not obtained
for lower Deborah number in Fig. 6 (b) and (d), where the
elastic effects are negligible, and no encapsulation was
obtained in purely viscous simulation of the flow in the
literature [3, 4], the encapsulation phenomena is also not
caused solely by the difference in the viscosity of the two
polymers.
At this point, we believe that polymer encapsulation in
coextrusion is a surface phenomenon, which is caused by
the difference in the wettability (spreading tendency of the
polymer on die surface) and surface tension of the two
polymers used. The polymer with higher wettability is
expected to encapsulate the other polymer during
coextrusion.
Conclusions
The Giesekus model was used for viscoelastic
simulation of a bi-layer flow in a square coextrusion die.
In contrast to the findings reported in the literature, in the
present work even with viscoelastic effects included in the
simulation, encapsulation of a high viscosity polymer by a
lower viscosity polymer, which is often observed in
coextrusion experiments, could not be captured. Since the
viscous encapsulation could not be captured with a purely
viscous formulation in our earlier work, and similar work
by other groups in the literature, it is concluded that
encapsulation is not caused by viscous or viscoelastic
effects. Instead, it is concluded that difference in
wettability and surface tension of the two polymers is
probably the major factor contributing to the
encapsulation phenomenon.
References
1. C. Rauwendaal, “Polymer Extrusion”, 4th Edition,
Hanser Publishers, (2001).
2. J. Dooley and L. Rudolph, J. of Plastic Film and
Sheeting, Vol. 19, 111 – 122 (2003).
3. A. Karagiannis, A. N. Hyrmak, and J. Vlachopoulos,
Rheologica Acta, Vol. 29, 71 – 87 (1990).
4. M. Gupta, SPE ANTEC Technical Papers, Vol. 54,
217 – 222 (2008).
5. K. B. Sunwoo, S. J. Park, S. J. Lee, K. H. Ahn, and S.
J. Lee, Rheologica Acta, Vol. 41, 144 – 153 (2002).
6. R. G. Larson, “Constitutive Equations for Polymeric
Melts and Solutions”, Butterworths, Boston, MA,
(1988).
7. J. N. Reddy, “An Introduction to Finite Element
Method”, 3rd Ed., McGraw Hill, (2006).
8. A. N. Brooks and T. J. R. Hughes, Computer
Methods in Applied Mechanics and Engineering,
Vol. 32, 199 – 259 (1982).
9. R. Guienette and M. Fortin, J. of Non-Newtonian
Fluid Mechanics, Vol. 60, 27 – 52 (1995).
10. T. Coupez and S. Marie, International Journal for
Supercomputer Applications and High Performance
Computing, Vol. 11, 277 – 285 (1997).
11. M. Gupta, SPE ANTEC Technical Papers, Vol. 56,
2032 – 2036 (2010).
Acknowledgement
This work was supported by the grant # 1046495 from the
National Science Foundation.
Fig. 1 Viscosity of the two polystyrenes.
Fig. 2 Interface shape in the experiments. The distance from the square entrance is (a) 0.4 L, (b) 2.42 L, and (c) 6.71 L, with L being the length of the side of the square [3].
Fig. 3 Interface shape predicted by a purely viscous simulation. The distance from the square entrance is (a) 0, (b) 0.42 L, and (c) 2.42 L, with L being the length of the side of the square [3].
Fig. 4 Axial velocity distribution in a bi-layer flow in a square channel for De=1.
Fig. 5 Secondary flow vortices in a bi-layer flow in a square channel for De=1.
Fig. 6 Shape on the interface in a bi-layer coextrusion in a square channel. Interface shape for (a) De=1, (b) De=0.1. Interface shape at channel exit for (c) De=1, (d) De=0.1.
Fig. 7 Magnitude of the second normal stress difference in a bi-layer flow in a square channel for De=1.
Fig. 8 Pressure distribution in a bi-layer flow in a square channel for De=1.
Return to
Paper of the Month.