# Numerical Simulation of Pressure and Velocity Profiles in Kneading Elements of a Co-Rotating Twin Screw Extruder.

V.L. BRAVO [1]

A.N. HRYMAK [*]

The objective of this work is to validate, via comparison with available experimental data, the results obtained from the numerical simulation of polymer melt flow in the kneading disc section of an intermeshing co-rotating twin screw extruder. A quasi-steady state 3-D solution of the conservation equations via the finite element method was obtained, and comparisons were made with experimental pressure profiles measured by McCullough and Hilton [1] on various kneading block elements. These measurements helped provide understanding of the flow patterns developed within the unit and provided a comprehensive approach of validating the numerical model. Results confirm the importance of a fully 3-D model for this type of geometry, where the model predicts the development of flow patterns in the radial directions and within the intermeshing region. The influence of inlet and outlet boundary conditions was studied and it was determined that they play an important role in the physical significance of the model solut ion. Comparisons of the simulation results with experimental data by McCullough and Hilton [1] for two different configurations of kneading discs showed good agreement, with some differences in the peaks of pressure produced at the narrow clearances encountered in intermeshing co-rotating twin screw extruders. Differences between simulation and experiments are attributed to a number of factors. It is difficult to measure the very steep pressure gradients generated over small lengths. The assumptions of isothermal flow and quasi-steady state may cause an over-prediction of the pressure peaks. Simulation results describe the general trends and produce good quantitative agreement in most of the kneading disc region.

1. INTRODUCTION

1.1. Intermeshing Co-rotating Twin Screw Extruder (ICRTSE)

Twin screw extruders include many machine configurations with widely different characteristics. In general, a twin screw extruder can be described as two parallel screws placed in a figure-eight-shaped barrel section. The relative position of the screws and their direction of rotation define different types of TSE. Among all the different types of twin screw extruders, the intermeshing co-rotating configuration is of special importance owing to its high efficiency in mixing, elimination of stagnant zones, and self-wiping action of the elements. The study of flow behavior within the complex geometry of a twin screw extruder is carried out using a numerical model. The optimization of the design parameters in any mixing system requires an understanding of the flow patterns developed within the unit. Using simplified mathematical models or experimental data and information obtained from extruders of similar design is a common practice in calculation procedures and parametric analyses for the improved design. In the particular case of the kneading disc section of a twin screw extruder, the intricate geometry, and the continuous motion of the mixing elements make the formulation and solution of mathematical models extremely complicated. Simplified models give approximate solutions that, in some cases, are not indicative of the true system response.

1.2 Mathematical Models and Experimental Work

From the theoretical model developed by Kaplan and Tadmor [2] for non-intermeshing machines based on the "three parallel plate" model to the modem 3D computational fluid dynamics (CFD) models, there has been a great deal of work done in the area of simulation of flow in twin screw extruders. The driving force for all this work has always been the understanding of the flow phenomena and the development of a tool for the prediction of mixing performance. Analytical models are limited to the prediction of pressure drop in regular conveying elements based on given flow rates. The use of geometrical correction factors is necessary to incorporate the effect of one screw on the drag and pressure terms of the other screw.

A more rigorous mathematical description of the complex flow passages encountered in the extruders is necessary when kneading discs are simulated. Tadmor et al. [3] used a control volume method (flow analysis network, FAN) for solving 2.5 dimensional flow problems in geometries with narrow gaps of variable thickness and different shapes. The FAN technique was used later by Szydlowski et al. [4] on the unwrapped flow channel of an ICRTSE neglecting variations in the radial direction with the use of the lubrication approximation.

Booy [5] described the cross section geometry of fully wiped co-rotating twin screws. In his work, Booy calculated the free volume between screws and barrel surface areas. The cross section of a twin screw has a unique shape for a given diameter, centerline distance, and number of parallel channels because of the requirement that one screw must wipe its mate in any position. Booy [6] described a mathematical flow model, with the assumption of stationary screw and rotating barrel, using the hypothetical case of a flexible barrel that moves around the screws with a velocity [pi][D.sub.S]N following the contour of the figure eight-shaped barrel.

Considerable work has been done in the 3-D simulation of flow in the kneading disc region of intermeshing co-rotating twin screw extruders. Gotsis et al. [7] proposed a model for the isothermal flow of a Newtonian fluid with minimal assumptions to geometry, using a finite element discretization. The objective is to simulate realistic flow patterns in the radial and axial directions. Mathematical complications arise as a consequence of considering the actual geometry. The flow pattern is time dependent because the flow boundaries change as the kneading discs rotate. To solve the problem of the time dependence of the geometry of the flow, they constructed a 3-D mesh covering the whole space inside the barrel that may be occupied either by the fluid or by the kneading discs. At each time step, the elements of the mesh are checked to determine whether they correspond to volumes occupied by the fluid or volumes occupied by the kneading blocks. The part of the space that is found to correspond to the fluid is assu med to obey the conservation equations whereas the part that is occupied by the kneading blocks during that time step is forced to rotate as a rigid body [7]. Lawal and Kalyon [8] used the same procedure as Gotsis et al. [7] to solve the conservation equations and extended the work for the analysis of intensity of segregation using particle tracers. Wang et al.[9] made calculations for the flow of viscous power law fluids in the various modules of a Krupp Werner & Pfleiderer ZSK-30 extruder. Kiani et at [10, 11] used the finite element method (FEM) and the spectral element method [12] for the 3-D quasisteady solution of the conservation equations for mass and momentum using periodic boundary conditions for the ZSK-30. Yang [13] developed 3-D simulations of the flow in the kneading disc sections of a ZSK-30 using the FEM and applied normal stress difference boundary conditions. The transient nature of the problem was addressed by Manas-Zloczower and Cheng [14]. Steady state flow calculations are made for each time step, setting the boundary conditions according to the no-slip condition at the walls and calculating the velocities at the surfaces of the kneading blocks based on the radius and the shaft rotation speed. The assumption of creeping flow is based on the fact that polymer melts have very high viscosities and very low Reynolds numbers (of order [10.sup-3] or less). The term Du/Dt is neglected from the momentum equations and the transient flow can be reasonably approximated as a sequence of steady state flows at intermediate times. Under this situation, the velocity and pressure fields adjust instantaneously relative to the rate at which the boundary configuration changes and consequently always appear to be in steady state with respect to the current configuration. However, the geometry configuration still depends on the change of the boundary position with respect to time. Yang [13] used the same approximation, also called the quasi-steady state solution. The boundary conditions imposed were the no-slip c ondition for the chamber wall and rotor surfaces, and a nominal value for the normal stress difference at the entrance and exit planes. Kajiwara et al. [15] simulated the flow in full flight elements using 3D finite elements, and carried out comparisons between co-rotating and counter-rotating twins screw extruders based on velocity, stresses and flow rates.

Rauwendaal [16] conducted studies on residence time distribution and mixing characteristics in ICRTSE and concluded that the overall extruder performance is dominated by the effect of the intermeshing region. Studies of twin screw extruders must include the flow behavior and mixing characteristics of the intermeshing region. This conclusion is widely accepted and has been the main motivation for moving from 2-D models into more rigorous 3-D models.

One of the major challenges in the study of twin screw extruders is the experimental validation of the results obtained by numerical models. McCullough and Hilton [1] presented an experimental investigation of the performance of co-rotating twin screw extruders. The data provided by McCullough and Hilton is not limited to measurements of pressure and temperature at the die, but include measurements along the barrel at different axial and radial positions. Christiano and Lindenfelzer [17] presented measurements of dynamic pressure distributions in kneading discs, with the feature of software developed to "remap" the radial pressure data files into a single file containing the entire pressure distribution over the length of the mixing element for a given screw position. These experimental results can be used in the future for further validation of the simulation results.

2. MODELING THE FLOW BEHAVIOR IN THE KNEADING DISC SECTION

2.1 The Kneading Disc Configuration

The intermeshing co-rotating twin screw extruder is one of the most extensively used pieces of equipment for the compounding and blending of polymers. The screws are made up of individual kneading blocks and screw elements that are assembled on the shafts to create the shear and mixing action required for a given material. The screw profile is designed so that the tip of one screw wipes the flank and root of the other screw, resulting in a self-wiping action.

The availability of information for this particular machine allows us to make comparisons and facilitates the verification of the validity of our results.

2.2 Geometry

The generation of the complex geometry involved in the kneading discs section of a twin screw extruder is based on the work done by Booy [5]. The calculation of the flow field within a twin screw extruder requires a complete knowledge of the geometry. The cross section of a twin screw has a unique shape given by diameter, centerline distance, and number of parallel channels because of the requirement that one screw must, in any position, wipe its mate.

The dimensions used for the computer geometry model are based on those of the Werner & Pleiderer ZSK-30 type co-rotating extruder, shown in Table 1. A stagger angle ([phi]) of 45[degrees] in a forward configuration was used for the simulations. Figure 1 describes the geometry of the cross section of kneading discs. The shape of the surface of one disc was generated using the equations developed by Booy [5]. The coordinates describing the surface of the disc were subsequently rotated by the stagger angle and translated to produce the rest of the discs. The resulting 3-D geometry represents the free volume between the kneading discs and the twin barrel. The geometry of the kneading discs is shown in Fig. 2. The transient nature of the flow problem in the kneading discs of an ICRTSE requires the re-generation of the geometry and the finite element mesh for every position of the rotors. Gotsis et al. [7] proposed a method where a single three-dimensional mesh covered the entire volume inside the barrel that may be occupi ed either by the fluid or by the kneading discs. This approach allowed Gotsis et al. [7] to eliminate the need for re-generating the finite element mesh for every position of the discs thus reducing one of the main complications in the solution of the time-dependent problem. The approach, although very appealing, has its limitations in the inaccuracies of the description of the rotor surfaces and especially the narrow gaps encountered in these systems.

The re-generation of the finite element mesh can be a tedious and time-consuming process. In this work, an automatic mesh generator was developed in order to overcome the problems associated with the multiple re-generation. Seven sequential geometries were specified to represent one quarter of the complete cycle of rotation for a geometry of five kneading discs in a forward configuration with a 45[degrees] stagger angle, according to Fig. 3. Because of symmetry, only one quarter of the complete cycle is necessary to represent the complete sequence of geometries. Variable [alpha] represents the angle between the left rotor tip and the y-axis for the first disc. Increments of 12.86[degrees] for each time step were taken, corresponding to a time Interval of 2.1429 x [10.sup.-2] seconds for a rotation speed of 100 rpm.

2.3 Governing Equations

Flow within the kneading discs region in a co-rotating twin screw extruder is modeled assuming the flow is in steady state, isothermal, creeping, incompressible and with no body forces.

The governing equations are the momentum equations without inertia terms and the continuity equation, which can be written, in dimensionless tensorial form as:

[delta]*[tau] - [delta]P = 0 (1)

[delta]*u = 0 (2)

where [tau] is the extra stress tensor, P the pressure and u the velocity vector. A Generalized Newtonian constitutive equation (18) is used to model the stress-strain relation for polymers:

[[tau].sub.ij] = [eta]([delta][u.sub.i]/[delta][x.sub.j] + [delta][u.sub.j]/[delta][x.sub.i]) (3)

The viscosity dependence on shearing is fitted with a Carreau model as:

[eta]([gamma]) = [[micro].sub.0][[1 + [([lambda][gamma]).sub.2]].sub.(n-1)/2] (4)

with

[gamma] = [square root of]1/2 [[pi].sub.[gamma]] (5)

[[pi].sub.[gamma]] = [[gamma].sub.ij] [[gamma].sub.ji] (6)

[[gamma].sub.ij] = [delta][u.sub.j]/[delta][x.sub.i] + [delta][u.sub.i]/[delta][x.sub.j] (7)

where [[micro].sub.0],[lambda] and n are rheological property parameters. The dimensionless groups used in Eqs 1 to 7 are:

x = x'/D u = u'/U

P = P'D/[[micro].sub.0]U [tau] = [tau]'D/[[micro].sub.0]U

[micro] = [micro]'/[[micro].sub.0] (8)

where D is the characteristic diameter (in this case, the outer diameter of the screw), U is the total average velocity (total flow rate over total area), [[micro].sub.0] is the zero-shear viscosity of the fluid and the prime (') denotes a dimensional variable.

2.4 Boundary Conditions

In a flow problem, the equations of conservation of mass and momentum must be solved with the appropriate boundary conditions in order to obtain meaningful results. At the barrel wall and the rotor surfaces, the standard no-slip boundary condition is used. The fluid at the barrel surface remains stationary whereas the fluid at the rotor surfaces moves according to the local rotation speed. The surface linear velocity can be calculated as the product of the angular velocity times the radius between the nodal point and the rotation axis as:

u(x, y, z, t) = [omega] X r(x, y, z, t) (9)

where [omega] is the angular velocity of the rotors and r is the position vector with respect to the center of each rotor.

The boundary conditions at the disc section inlet and outlet are more complicated to address because the section is being isolated from the upstream and downstream sections. As Yang [13] points out, it is very difficult to calculate a priori the true velocity profiles at the entrance and exit planes. The kneading disc sections are usually fed by screw conveying sections. The boundary condition at the inlet of the kneading disc section should be the exit values of the preceding conveying section, and it is expected that both the feeding and the entered sections affect the velocity profiles at the meeting planes. The complete specification of the inlet velocity profile would require the calculation of the flow field in the screw conveying section in conjunction with the kneading disc section. Kianni and Samman [11] treated the boundary velocity condition between the two axial surfaces as periodic. Gotsis et al. [7] imposed pressure conditions at the entry and exit of the kneading disc section. Yang [13], used a nominal value for the normal stress difference, combining pressure and viscous stresses, in the axial direction between the two planes.

Bravo [19] presented an alternative approach to the inlet/outlet boundary condition by imposing fully developed velocity profiles at the inlet and outlet at a specified distance from the rotating discs. This approach is limited to conditions where the maximum pressure occurs at the entrance and a pressure drop exists throughout the kneading disc section [1]. The normal stress boundary condition, combining pressure and viscous stress, is used in this work. The main advantages of the normal stress approach are: a) the ability to impose positive and negative pressure gradients between inlet and outlet, allowing the calculation of flow rate as a function of pressure and drag flows; and b) the elimination of the need for calculating inlet and outlet velocities a priori.

Application of the Galerkin method to the momentum equation produces:

[[R.sup.i].sub.m] = [[integral of].sub.[omega]] {[nabla]*[sigma]}[N.sup.i] d[omega] = [[integral of].sub.[omega]] {[nabla]*[N.sup.i]}[sigma]d[omega] - [[integral of].sub.[gamma]] n*[sigma][N.sup.i]d[gamma] (10)

where the first term on the right hand side applies to the entire domain, while the second term applies to the boundaries. Let us consider the boundary term,

[F.sub.i] = [[integral of].sub.[gamma] n*[sigma][N.sup.i] d[gamma] (11)

This term can be expanded to obtain:

[F.sub.x] = [[integral of].sub.[gamma]] [n.sub.x] [[sigma].sub.xx][N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.y] [[sigma].sub.xy] [N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.z] [[sigma].sub.xz] [N.sup.i] d[gamma] (12)

[F.sub.y] = [[integral of].sub.[gamma]] [n.sub.x] [[sigma].sub.yx][N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.y] [[sigma].sub.yy] [N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.z] [[sigma].sub.yz] [N.sup.i] d[gamma] (13)

[F.sub.z] = [[integral of].sub.[gamma]] [n.sub.x] [[sigma].sub.zx][N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.y] [[sigma].sub.zy] [N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.z] [[sigma].sub.zz] [N.sup.i] d[gamma] (14)

In the above equations, [F.sub.i] represents the resultant forces and [N.sup.i] represents the shape functions. Also [N.sub.1], [n.sub.2] and [n.sub.3] are the components of the normal unit vector at the boundary, which for plane [eta]-[zeta] in the normal coordinate system [24] are:

n = [[a.sub.1], [a.sub.2], [a.sub.3]]/[square root of][[a.sup.2].sub.1] + [[a.sup.2].sub.2] + [[a.sup.2].sub.3] (15)

with

[a.sub.1] = [delta]y/[delta][eta] [delta]z/[delta][zeta] - [delta]y/[delta][zeta] [delta]z/[delta][eta] (16)

[a.sub.2] = [delta]x/[delta][zeta] [delta]z/[delta][eta] - [delta]x/[delta][eta] [delta]z/[delta][zeta] (17)

[a.sub.3] = [delta]x/[delta][eta] [delta]y/[delta][zeta] - [delta]x/[delta][zeta] [delta]y/[delta][eta] (18)

A normal stress difference between the entrance and exit surfaces was set as [delta] [[sigma].sub.xx], such that the desired flow rate was obtained. Equations 12, 13 and 14 are applied over the inlet and outlet boundaries. The average pressure rise measured at the apex ports over the length of 42 mm in the 45/5/42 configuration was approximately 650 kPa, and the average pressure rise measured in the side ports was 500 kPa. The reference value of [[sigma].sub.xx] was set zero at the entrance plane and [delta][[sigma].sub.xx] was set to 1035 kPa to match the experimental flow rate of 4.35 Kg/h. the rotation speed was either 60 or 140 RPM.

The problem of the formulation of appropriate boundary conditions at inflow and outflow planes has been a controversial one. Sometimes it is possible to obtain seemingly acceptable solutions for the finite flow field that are, in reality, poor approximations of the real problem [20]. Validation via experimental flow visualization and velocity measurements is the best way to assess model boundary conditions.

2.5 Finite Element Formulation

A three-dimensional finite element method is used to discretize the governing Eqs 1 and 2 in the kneading disc region. Several texts provide the necessary mathematical background of the method applied to fluid mechanics problems [21, 22]. The field variables [velocities and pressures] are approximated by 27-node triquadratic brick elements for velocities and 8 node trilinear brick elements for pressure. For the discretization of the twin screw domain, different densities of 27-node triquadratic brick element meshes were tested. The results shown in this work correspond to a mesh of 2,016 elements and 20,881 nodes, as depicted in Fig. 4.

The finite element code was used to obtain the velocity profiles and the pressure distributions of the flow field according to experimental stagger angles and operating conditions by McCullough and Hilton [1].

3. COMPARISON WITH EXPERIMENTAL PRESSURE PROFILES

Traditional twin screw extrusion process data are limited to pressures and temperatures measured at the die. Measurements are not normally made along the barrel since the best measurement locations are a function of screw configuration. Steep pressure gradients are likely to appear in both the axial and rotation directions. McCullough and Hilton [1] developed an experimental apparatus consisting of a 30 mm diameter TSE using a 2L/D feed adapter, an instrumented 6L/D test barrel section and standard die assembly. The configuration of the apparatus is shown in Fig. 5. The reverse kneading elements at the end of the kneading disc section restricted the flow and raised the pressure within the system. Figure 5a shows an end view of the barrel section from the feed end with the location of the pressure transducer at the top of the apex and the side of the main screw on the right. To obtain the required pressure data, McCullough and Hilton equipped the barrel with five ports evenly spaced every 30 mm along the top and bottom of the center apex region, and three ports along the sides spaced every 60 mm. The TSE assembly was mounted on a movable lathe bed. This allowed the barrel section to be moved back and forth along the screw elements for the various measurements. High-speed data measurements confirmed that the pressure transducers were capable of responding as quickly as 800 MPa/s.

Response time is critical because at 100 rpm (or 600 degrees/s) pressure can fluctuate as much as 8 MPa in 0.05 seconds, or 160 MPa/s. The data collection frequency was set to record pressure data at each degree of screw rotation: 360 records/second at 60 rpm and 800 records per second at 140 rpm.

Results are presented for runs at 60 and 140 rpm using a low density polyethylene. Pressures were measured at the side of the main screw and at the top of the apex, for different arrays of kneading discs, according to Fig. 5b. Of special interest for the comparisons with simulation data are two of the configurations depicted in Fig. 5b. The configuration of forward kneading discs are 45/5/20 (where 45 is the stagger angle, 5 is the number of discs and 20 stands for a total disc section length of 20 mm) and 45/5/42. The pressure measurements are plotted as a function of the amount of rotation. The axial pressure drop plots show average, maximum and minimum values over one cycle of rotation for apex and side measurements. For the simulations, the average side pressure drop was used as an initial guess for the normal stress difference imposed as boundary condition and then adjusted to obtain the correct flow rate. The average side pressure was used since it is the most representative of the average behavior in the system. The pressure fluctuations in the apex region are more intense than at the side because of the interaction of the two screws at the intermeshing region. The flow rate was kept at 4.35 kg/hr to match experimental conditions The complete set of operating conditions and material properties (1) are presented in Tables 2 and 3 and Fig. 6.

A total of 7 steady state solutions with increments of 12.86[degrees] are used to complete one quarter of a revolution. Because of the symmetry of the system, calculation of only one quarter of a revolution is necessary to describe the complete cycle. The complete set of solutions is obtained by means of spatial transformation to rotate the first quarter geometry 180[degrees] around the x-axis to obtain the second quarter. Once the solution for a half revolution is obtained, symmetry applies again making the second half revolution identical to the first half. An isometric view of the combination of kneading blocks used for the experimental work is shown in Fig. 7.

The numerical determination of pressures at the side and apex is done by interpolating the values at the locations equivalent to the location of the pressure transducers in the experimental setup. The reported experimental pressure value is centered at a point, measured over a small area. Nine calculated pressure values are interpolated in a cross pattern as shown in Fig. 8 to cover a 1-mm-diameter circle, which corresponds to the area exposed to the pressure transducer in the experimental setup. An average of the nine points is then taken to account for the pressure on the 1-mm-diameter circle. The pressure values are interpolated from the nodal values previously obtained through the FEM calculations. Once the coordinate of the point where the pressure needs to be calculated is determined, an algorithm searches for the element of the mesh that contains that coordinate. The pressure is then interpolated from the nodal values using the element's interpolation functions. For a full revolution, the numerical pr ocedure of determining the pressure at the apex and side ports is repeated 28 times for every disc.

For the simulation of the 45/5/42 configuration the geometry depicted in Fig. 5 representing 5 pairs of discs was used. For the 45/5/20 configuration a system with two pairs of discs was used. Two disks were used for the 45/5/20 configuration to save computer time, given the limited computational resources. However, from Fig. 5b, there is an important effect of disk 4 of the 45/5/20 block on the circumferential pressure profile of disk 5, due to their proximity. Disk 1 of the 45/5/42 is aligned with disk 5 of the 45/5/20 block, and disk 2 of the 45/5/42 block does not have a major effect. Experimental and simulation results have a better match for this particular case, because of the distance of disk 2 of 45/5/42 from disk 5 of 45/5/20. However, a more complete analysis would include more disks in the 45/5/20 configuration.

4. DISCUSSION

Figures 9 and 10 show the comparison between simulation and experimental variation in pressure for a complete cycle of disc 5 in the 45/5/20 configuration. At this location, the pressure fluctuations are the largest in the system. The experimental peaks for the apex port reach 4250 kPa whereas the minimum pressure is 1750 kPa. Figure 10 shows the pressure profiles at the side port, where the amplitude of the fluctuations is less than for the apex port. Figures 11 to 13 present experimental and simulation results for the individual discs of the 45/5/42 configuration at the apex port. Figures 14 to 16 show the experimental and simulation results at the side port of the 45/5/42 configuration.

An increasing pressure in the axial direction has been reported [1] for the configuration of forward kneading blocks, as predicted by the simulations. Although the staggering of the discs reduces the conveying capacity of kneading blocks, they still maintain some ability to move material forward. The forward conveying is produced by the dragging effect whereas the axial pressure induces a back flow component that ensures a fully filled channel and hence better mixing performance. The normal stress boundary condition determines whether there is a positive or negative conveying effect for specific pressure drop and rotation speed conditions.

The motion of the rotors causes the fluctuating pressures observed experimentally and predicted with the model. Different mechanisms cause the pressure fluctuations in the side and apex. At the side observation point, the advancing tip of the discs creates a converging region between the disc and the barrel causing high-pressure peaks. The relative motion of the trailing side of the disc to the barrel wall produces the low pressures, as illustrated in Fig. 17. The high pressures at the apex originate from the encounter of the tips of the two discs. A small chamber is created and then reduced in size until it completely disappears. squeezing the fluid and creating the high peaks of pressure. Immediately after the chamber disappears, a small chamber appears again and increases in size until it communicates with the wide channels at both sides of the chamber, creating the low pressure values, as illustrated in Fig. 18.

Pressure valleys at one disc coincide with the highest pressure peaks in the subsequent disc, generating the experimentally observed back flows at the rear side of the discs. This effect is very much affected by the stagger angle as well as the thickness of the discs. It is expected that thicker discs would generate wider valleys of pressure and therefore greater back flows. The 45[degrees] stagger angle kneading blocks have the particular property of the coincidence of the highest pressure peaks on the advancing side of the disc's tip with the lowest pressure valleys on the rear side of the contiguous disc's tip. Based on this observation and from the model calculations, the phenomena observed by McCullough and Hilton [1] can be explained. Figures 9 and 10 show the pressure behavior as a function of the angle of rotation for disc 5 in the 45/5/20 configuration. Disc 5 is adjacent and aligned with the first disc of the 45/5/42 configuration (Fig. 5b). The resulting effect is that of a long disc by adding the 4 mm thickness of the 45/5/20 disc 5 to the 8.2 mm of the 45/5/42 disc 1, for 12.2 mm total thickness. Simulations reveal that the pressure peaks for a single pair of discs, under similar flow rate and rotational speed to those used in the experimental work, are higher than the peaks produced in discs preceded and followed by another disc as occurs in real kneading blocks. This behavior has its explanation in the fact that the peaks and valleys of pressure for one disc are dampened by the effect of discs behind and ahead of it. In the case of disc 5 of the 45/5/20, the subsequent downstream staggered disc is at a distance of 8.22 mm, thus reducing its effect over the pressure profile of the preceding disc. This explains the experimental observation of an increase in the pressure fluctuations from the peaks to the valleys when approaching disc 5. It also explains why the amount of fluctuation decreases rapidly once the 45/5/42 configuration is reached. The intensity of pressure fluctuations is affected by the reverse elements located at the end of the 45/5/42 configuration that actually holds up the polymer melt and creates an almost uniform pressure distribution over the cross section of disc 5 of the 45/5/42 block. The simulation of the 45/5/42 configuration (Fig. 5b) used 5 pairs of discs. For the 45/5/20 configuration, two pairs of discs were used in the simulation. Disc 5 of the 45/5/20 block is aligned with the first disc of the 45/5/42 block, thus its local region is different from the rest of the discs in the same configuration. Comparisons of the model with experimental data on the 45/5/20 configuration (Figs. 9 and 10) show very good correlation, and together with the simulations for the 45/5/42 configuration help to explain the pressure behavior observed experimentally. Simulations for the 45/5/42 show also good correlation for all discs but disc 5, which is expected since it is adjacent to the reverse kneading blocks.

Although the model predictions are reasonably close to experimental pressure data, there is a consistent difference in the pressure peaks. When analyzing the model, there are two main assumptions made among those discussed in Sections 2 and 3: quasi-steady state and isothermal behavior. The quasi steady state assumption is said to apply for systems where Re (Reynold number, Re = [rho]VD/[micro]) [less than][less than]1 and Sr (Strouhal number Sr = [omega]D/V) [less than][less than]1. The assumption of Re [less than][less than]1 for polymeric flows is valid due to the very high viscosity of polymer melts and hence acceleration terms are neglected. The assumption of Sr [less than][less than]1 implies that the time scale of the perturbation, i.e. a periodic change in boundary conditions such as a pulsating flow or a moving boundary, must be much larger than the time scale of the diffusion of momentum. i.e. t[greater than][greater than][D.sup.2]/[nu]. This assumption is valid for polymeric flows that demonstrate ve ry low values of [D.sup.2]/[nu]. There are very large pressure gradients generated in a very small time interval due to the sudden transition from the narrow gap of the disc-barrel clearance to the wider channel that follows. These pressure peaks represent extremely small zones within the flow field that may not affect flow quantities for purposes such as determining mixing performance from numerically calculated flow fields. Inclusion of transient terms and an algorithm to integrate the equations in time may capture the sharp pressure transitions accurately.

The system is assumed to be isothermal. In real systems, it is expected that zones of high shear would experience higher temperatures due to viscous dissipation and therefore lower values of viscosity than estimated from the isothermal model. It is well known that the zones where the pressure peaks occur are also the zones with the highest shear rates in the system. The isothermal flow assumption could be a reason for overestimation of the pressure peaks in the system.

An additional reason for the discrepancies between experimental and simulation pressures could be the model for the dependence of viscosity on the magnitude of the rate-of-strain tensor [gamma] (see Eqs 4 to 7). The highest peaks of pressure occur at the narrow gaps where the [gamma] values are the highest. The Carreau model fits the data for a range of [gamma] between 0 and 5000 [s.sup.-1] according to Fig. 6. Calculations have demonstrated that at the narrow gaps the values of [gamma] can be greater that 5000 [s.sup.-1]. Over-prediction of the viscosity at high values of [gamma] could produce an overestimation of the pressure profiles. An alternative to the present model, would be a piece-wise fit of the Carreau model to specified ranges of shear rate to capture the apparently lower power law index at shear rates in the 1000-10,000 [s.sup.-1] range.

The results are also in qualitative agreement with circumferential pressure measurements presented by Kiani and Heidemeyer [23] and calculations made by Szydlowski et al. [4]. Physical limitations in the pressure transducers--i.e., speed of response, response to a non-uniform pressure in the face of the transducer, noise, etc.--may also introduce deviations from the prediction made by the model. Pressure transducers have limitations in the measurement of steep gradients over small areas since they are usually designed to measure constant pressures over the area of the probe. Limitations in the speed of response of the transducer could produce a "dampening" of the signal at those locations where the gradients are the steepest.

5. CONCLUSIONS

A quasi-steady state 3-D finite element model has been implemented and compared against experimental data provided by McCullough and Hilton [1]. For the purpose of comparison, simulation results were analyzed over small areas equivalent to that of the pressure transducers, since numerical interpolation determines the pressure value at a point. Differences are noticed at the pressure peaks produced by the tips of the kneading discs "sweeping" the surface of the barrel, whereas for the other regions the comparison is reasonably good. Sources of error to explain the differences between simulation and experiment include: mesh density, isothermal assumption, fit of the constitutive equation at high shear rates, and pressure transducer response time.

ACKNOWLEDGMENTS

The authors wish to acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada and Investigacion y Desarrollo C.A. (Venezuela).

(1.) V. L. Bravo on leave from Investigacion y Desarrollo C.A. (INDESCA), P.O. Box 10319. Complejo Petroquimico El Tablazo. Maracalbo, 4001. Venezuela.

(*.) Author to whom correspondence should be addressed

(email: hrymak@mcmaster.ca, fax: (905) 521-1350).

6. REFERENCES

(1.) T. W. McCullough and B. T. Hilton, SPE ANTEC Technical Papers, 3372 (1993).

(2.) A. Kaplan and Z. Tadmor, Polym. Eng. Sci., 14, 58 (1974).

(3.) Z. Tadmor, E. Broyer, and C. Gutfinger, Polym. Eng. Sci., 14, 660 (1974].

(4.) W. Szydlowski and J. L. White, Intern. Polym. Processing, 2, 142 (1987).

(5.) M. L. Booy, Polym. Eng. Sci., 18, 3 (1978).

(6.) M. L. Booy, Polym. Eng. Sci., 20, 1220 (1980).

(7.) A. D. Gotsis, Z. Ji, and D. M. Kaylon, SPE ANTEC Technical Papers. 139 (1990).

(8.) A. Lawal and D. M. Kalyon, Polym. Eng. Sci., 35, 1325 (1995).

(9.) Y. Wang. J. L. White, and W. Szydlowski, Intern. Polym. Proc., 4, 262 (1989)

(10.) A. Kiani, E. W. Grald, and S. Subbiah, "Spectral Element Analysis of the Mixing Characteristics of Twin Screw Extruders," Session on Numerical Simulation of Mixing Phenomena. AIChE Annual Meeting (1992).

(11.) A. Kiani and H. J. Samann, SPE ANTEC Technical Papers, 2578 (1993).

(12.) A. T. Patera, J. Computational Physics, 54, 468 (1984).

(13.) H. H. Yang, Flow Field Analysis of Batch and Continuous Mixing Equipment, PhD thesis. Case-Western Reserve University, Cleveland, Ohio (1993).

(14.) I. Manas-Zloczower and J. J. Cheng, Polym. Eng. Sci., 29, 1059 (1989).

(15.) T. Kajiwara, Y. Nagashima, Y. Nakano, and K. Funatsu, Polym. Eng. Sci., 36, 2142 (1996)

(16.) C. Rauwendaal, Polym. Eng. Sci., 21, 1092 (1981).

(17.) J. Christiano and M. Lindenfelzer, SPE ANTEC Technical Papers, 78 (1997).

(18.) R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena. John Wiley (1984).

(19.) V. Bravo, 3D Simulation of Quasi-Steady Flow Field in Twin Screw Extruders, M. Eng. thesis, McMaster University. Hamilton, Ontario, Canada (1995).

(20.) J. S. Vrentas, C. M. Vrentas, and H. C. Ling, J Non-Newt. Fluid Mech., 19, 1 (1985).

(21.) G. Dhatt and G. Touzot, The Finite Element Method Displayed, John Wiley and Sons, New York (1984).

(22.) J. N. Reddy, Introduction to the Finite Element Method, McGraw-Hill, New York (1984).

(23.) A. Kiani and P. Heidemeyer, SPE ANTEC Technical Papers, 94 (1997).

A.N. HRYMAK [*]

The objective of this work is to validate, via comparison with available experimental data, the results obtained from the numerical simulation of polymer melt flow in the kneading disc section of an intermeshing co-rotating twin screw extruder. A quasi-steady state 3-D solution of the conservation equations via the finite element method was obtained, and comparisons were made with experimental pressure profiles measured by McCullough and Hilton [1] on various kneading block elements. These measurements helped provide understanding of the flow patterns developed within the unit and provided a comprehensive approach of validating the numerical model. Results confirm the importance of a fully 3-D model for this type of geometry, where the model predicts the development of flow patterns in the radial directions and within the intermeshing region. The influence of inlet and outlet boundary conditions was studied and it was determined that they play an important role in the physical significance of the model solut ion. Comparisons of the simulation results with experimental data by McCullough and Hilton [1] for two different configurations of kneading discs showed good agreement, with some differences in the peaks of pressure produced at the narrow clearances encountered in intermeshing co-rotating twin screw extruders. Differences between simulation and experiments are attributed to a number of factors. It is difficult to measure the very steep pressure gradients generated over small lengths. The assumptions of isothermal flow and quasi-steady state may cause an over-prediction of the pressure peaks. Simulation results describe the general trends and produce good quantitative agreement in most of the kneading disc region.

1. INTRODUCTION

1.1. Intermeshing Co-rotating Twin Screw Extruder (ICRTSE)

Twin screw extruders include many machine configurations with widely different characteristics. In general, a twin screw extruder can be described as two parallel screws placed in a figure-eight-shaped barrel section. The relative position of the screws and their direction of rotation define different types of TSE. Among all the different types of twin screw extruders, the intermeshing co-rotating configuration is of special importance owing to its high efficiency in mixing, elimination of stagnant zones, and self-wiping action of the elements. The study of flow behavior within the complex geometry of a twin screw extruder is carried out using a numerical model. The optimization of the design parameters in any mixing system requires an understanding of the flow patterns developed within the unit. Using simplified mathematical models or experimental data and information obtained from extruders of similar design is a common practice in calculation procedures and parametric analyses for the improved design. In the particular case of the kneading disc section of a twin screw extruder, the intricate geometry, and the continuous motion of the mixing elements make the formulation and solution of mathematical models extremely complicated. Simplified models give approximate solutions that, in some cases, are not indicative of the true system response.

1.2 Mathematical Models and Experimental Work

From the theoretical model developed by Kaplan and Tadmor [2] for non-intermeshing machines based on the "three parallel plate" model to the modem 3D computational fluid dynamics (CFD) models, there has been a great deal of work done in the area of simulation of flow in twin screw extruders. The driving force for all this work has always been the understanding of the flow phenomena and the development of a tool for the prediction of mixing performance. Analytical models are limited to the prediction of pressure drop in regular conveying elements based on given flow rates. The use of geometrical correction factors is necessary to incorporate the effect of one screw on the drag and pressure terms of the other screw.

A more rigorous mathematical description of the complex flow passages encountered in the extruders is necessary when kneading discs are simulated. Tadmor et al. [3] used a control volume method (flow analysis network, FAN) for solving 2.5 dimensional flow problems in geometries with narrow gaps of variable thickness and different shapes. The FAN technique was used later by Szydlowski et al. [4] on the unwrapped flow channel of an ICRTSE neglecting variations in the radial direction with the use of the lubrication approximation.

Booy [5] described the cross section geometry of fully wiped co-rotating twin screws. In his work, Booy calculated the free volume between screws and barrel surface areas. The cross section of a twin screw has a unique shape for a given diameter, centerline distance, and number of parallel channels because of the requirement that one screw must wipe its mate in any position. Booy [6] described a mathematical flow model, with the assumption of stationary screw and rotating barrel, using the hypothetical case of a flexible barrel that moves around the screws with a velocity [pi][D.sub.S]N following the contour of the figure eight-shaped barrel.

Considerable work has been done in the 3-D simulation of flow in the kneading disc region of intermeshing co-rotating twin screw extruders. Gotsis et al. [7] proposed a model for the isothermal flow of a Newtonian fluid with minimal assumptions to geometry, using a finite element discretization. The objective is to simulate realistic flow patterns in the radial and axial directions. Mathematical complications arise as a consequence of considering the actual geometry. The flow pattern is time dependent because the flow boundaries change as the kneading discs rotate. To solve the problem of the time dependence of the geometry of the flow, they constructed a 3-D mesh covering the whole space inside the barrel that may be occupied either by the fluid or by the kneading discs. At each time step, the elements of the mesh are checked to determine whether they correspond to volumes occupied by the fluid or volumes occupied by the kneading blocks. The part of the space that is found to correspond to the fluid is assu med to obey the conservation equations whereas the part that is occupied by the kneading blocks during that time step is forced to rotate as a rigid body [7]. Lawal and Kalyon [8] used the same procedure as Gotsis et al. [7] to solve the conservation equations and extended the work for the analysis of intensity of segregation using particle tracers. Wang et al.[9] made calculations for the flow of viscous power law fluids in the various modules of a Krupp Werner & Pfleiderer ZSK-30 extruder. Kiani et at [10, 11] used the finite element method (FEM) and the spectral element method [12] for the 3-D quasisteady solution of the conservation equations for mass and momentum using periodic boundary conditions for the ZSK-30. Yang [13] developed 3-D simulations of the flow in the kneading disc sections of a ZSK-30 using the FEM and applied normal stress difference boundary conditions. The transient nature of the problem was addressed by Manas-Zloczower and Cheng [14]. Steady state flow calculations are made for each time step, setting the boundary conditions according to the no-slip condition at the walls and calculating the velocities at the surfaces of the kneading blocks based on the radius and the shaft rotation speed. The assumption of creeping flow is based on the fact that polymer melts have very high viscosities and very low Reynolds numbers (of order [10.sup-3] or less). The term Du/Dt is neglected from the momentum equations and the transient flow can be reasonably approximated as a sequence of steady state flows at intermediate times. Under this situation, the velocity and pressure fields adjust instantaneously relative to the rate at which the boundary configuration changes and consequently always appear to be in steady state with respect to the current configuration. However, the geometry configuration still depends on the change of the boundary position with respect to time. Yang [13] used the same approximation, also called the quasi-steady state solution. The boundary conditions imposed were the no-slip c ondition for the chamber wall and rotor surfaces, and a nominal value for the normal stress difference at the entrance and exit planes. Kajiwara et al. [15] simulated the flow in full flight elements using 3D finite elements, and carried out comparisons between co-rotating and counter-rotating twins screw extruders based on velocity, stresses and flow rates.

Rauwendaal [16] conducted studies on residence time distribution and mixing characteristics in ICRTSE and concluded that the overall extruder performance is dominated by the effect of the intermeshing region. Studies of twin screw extruders must include the flow behavior and mixing characteristics of the intermeshing region. This conclusion is widely accepted and has been the main motivation for moving from 2-D models into more rigorous 3-D models.

One of the major challenges in the study of twin screw extruders is the experimental validation of the results obtained by numerical models. McCullough and Hilton [1] presented an experimental investigation of the performance of co-rotating twin screw extruders. The data provided by McCullough and Hilton is not limited to measurements of pressure and temperature at the die, but include measurements along the barrel at different axial and radial positions. Christiano and Lindenfelzer [17] presented measurements of dynamic pressure distributions in kneading discs, with the feature of software developed to "remap" the radial pressure data files into a single file containing the entire pressure distribution over the length of the mixing element for a given screw position. These experimental results can be used in the future for further validation of the simulation results.

2. MODELING THE FLOW BEHAVIOR IN THE KNEADING DISC SECTION

2.1 The Kneading Disc Configuration

The intermeshing co-rotating twin screw extruder is one of the most extensively used pieces of equipment for the compounding and blending of polymers. The screws are made up of individual kneading blocks and screw elements that are assembled on the shafts to create the shear and mixing action required for a given material. The screw profile is designed so that the tip of one screw wipes the flank and root of the other screw, resulting in a self-wiping action.

The availability of information for this particular machine allows us to make comparisons and facilitates the verification of the validity of our results.

2.2 Geometry

The generation of the complex geometry involved in the kneading discs section of a twin screw extruder is based on the work done by Booy [5]. The calculation of the flow field within a twin screw extruder requires a complete knowledge of the geometry. The cross section of a twin screw has a unique shape given by diameter, centerline distance, and number of parallel channels because of the requirement that one screw must, in any position, wipe its mate.

The dimensions used for the computer geometry model are based on those of the Werner & Pleiderer ZSK-30 type co-rotating extruder, shown in Table 1. A stagger angle ([phi]) of 45[degrees] in a forward configuration was used for the simulations. Figure 1 describes the geometry of the cross section of kneading discs. The shape of the surface of one disc was generated using the equations developed by Booy [5]. The coordinates describing the surface of the disc were subsequently rotated by the stagger angle and translated to produce the rest of the discs. The resulting 3-D geometry represents the free volume between the kneading discs and the twin barrel. The geometry of the kneading discs is shown in Fig. 2. The transient nature of the flow problem in the kneading discs of an ICRTSE requires the re-generation of the geometry and the finite element mesh for every position of the rotors. Gotsis et al. [7] proposed a method where a single three-dimensional mesh covered the entire volume inside the barrel that may be occupi ed either by the fluid or by the kneading discs. This approach allowed Gotsis et al. [7] to eliminate the need for re-generating the finite element mesh for every position of the discs thus reducing one of the main complications in the solution of the time-dependent problem. The approach, although very appealing, has its limitations in the inaccuracies of the description of the rotor surfaces and especially the narrow gaps encountered in these systems.

The re-generation of the finite element mesh can be a tedious and time-consuming process. In this work, an automatic mesh generator was developed in order to overcome the problems associated with the multiple re-generation. Seven sequential geometries were specified to represent one quarter of the complete cycle of rotation for a geometry of five kneading discs in a forward configuration with a 45[degrees] stagger angle, according to Fig. 3. Because of symmetry, only one quarter of the complete cycle is necessary to represent the complete sequence of geometries. Variable [alpha] represents the angle between the left rotor tip and the y-axis for the first disc. Increments of 12.86[degrees] for each time step were taken, corresponding to a time Interval of 2.1429 x [10.sup.-2] seconds for a rotation speed of 100 rpm.

2.3 Governing Equations

Flow within the kneading discs region in a co-rotating twin screw extruder is modeled assuming the flow is in steady state, isothermal, creeping, incompressible and with no body forces.

The governing equations are the momentum equations without inertia terms and the continuity equation, which can be written, in dimensionless tensorial form as:

[delta]*[tau] - [delta]P = 0 (1)

[delta]*u = 0 (2)

where [tau] is the extra stress tensor, P the pressure and u the velocity vector. A Generalized Newtonian constitutive equation (18) is used to model the stress-strain relation for polymers:

[[tau].sub.ij] = [eta]([delta][u.sub.i]/[delta][x.sub.j] + [delta][u.sub.j]/[delta][x.sub.i]) (3)

The viscosity dependence on shearing is fitted with a Carreau model as:

[eta]([gamma]) = [[micro].sub.0][[1 + [([lambda][gamma]).sub.2]].sub.(n-1)/2] (4)

with

[gamma] = [square root of]1/2 [[pi].sub.[gamma]] (5)

[[pi].sub.[gamma]] = [[gamma].sub.ij] [[gamma].sub.ji] (6)

[[gamma].sub.ij] = [delta][u.sub.j]/[delta][x.sub.i] + [delta][u.sub.i]/[delta][x.sub.j] (7)

where [[micro].sub.0],[lambda] and n are rheological property parameters. The dimensionless groups used in Eqs 1 to 7 are:

x = x'/D u = u'/U

P = P'D/[[micro].sub.0]U [tau] = [tau]'D/[[micro].sub.0]U

[micro] = [micro]'/[[micro].sub.0] (8)

where D is the characteristic diameter (in this case, the outer diameter of the screw), U is the total average velocity (total flow rate over total area), [[micro].sub.0] is the zero-shear viscosity of the fluid and the prime (') denotes a dimensional variable.

2.4 Boundary Conditions

In a flow problem, the equations of conservation of mass and momentum must be solved with the appropriate boundary conditions in order to obtain meaningful results. At the barrel wall and the rotor surfaces, the standard no-slip boundary condition is used. The fluid at the barrel surface remains stationary whereas the fluid at the rotor surfaces moves according to the local rotation speed. The surface linear velocity can be calculated as the product of the angular velocity times the radius between the nodal point and the rotation axis as:

u(x, y, z, t) = [omega] X r(x, y, z, t) (9)

where [omega] is the angular velocity of the rotors and r is the position vector with respect to the center of each rotor.

The boundary conditions at the disc section inlet and outlet are more complicated to address because the section is being isolated from the upstream and downstream sections. As Yang [13] points out, it is very difficult to calculate a priori the true velocity profiles at the entrance and exit planes. The kneading disc sections are usually fed by screw conveying sections. The boundary condition at the inlet of the kneading disc section should be the exit values of the preceding conveying section, and it is expected that both the feeding and the entered sections affect the velocity profiles at the meeting planes. The complete specification of the inlet velocity profile would require the calculation of the flow field in the screw conveying section in conjunction with the kneading disc section. Kianni and Samman [11] treated the boundary velocity condition between the two axial surfaces as periodic. Gotsis et al. [7] imposed pressure conditions at the entry and exit of the kneading disc section. Yang [13], used a nominal value for the normal stress difference, combining pressure and viscous stresses, in the axial direction between the two planes.

Bravo [19] presented an alternative approach to the inlet/outlet boundary condition by imposing fully developed velocity profiles at the inlet and outlet at a specified distance from the rotating discs. This approach is limited to conditions where the maximum pressure occurs at the entrance and a pressure drop exists throughout the kneading disc section [1]. The normal stress boundary condition, combining pressure and viscous stress, is used in this work. The main advantages of the normal stress approach are: a) the ability to impose positive and negative pressure gradients between inlet and outlet, allowing the calculation of flow rate as a function of pressure and drag flows; and b) the elimination of the need for calculating inlet and outlet velocities a priori.

Application of the Galerkin method to the momentum equation produces:

[[R.sup.i].sub.m] = [[integral of].sub.[omega]] {[nabla]*[sigma]}[N.sup.i] d[omega] = [[integral of].sub.[omega]] {[nabla]*[N.sup.i]}[sigma]d[omega] - [[integral of].sub.[gamma]] n*[sigma][N.sup.i]d[gamma] (10)

where the first term on the right hand side applies to the entire domain, while the second term applies to the boundaries. Let us consider the boundary term,

[F.sub.i] = [[integral of].sub.[gamma] n*[sigma][N.sup.i] d[gamma] (11)

This term can be expanded to obtain:

[F.sub.x] = [[integral of].sub.[gamma]] [n.sub.x] [[sigma].sub.xx][N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.y] [[sigma].sub.xy] [N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.z] [[sigma].sub.xz] [N.sup.i] d[gamma] (12)

[F.sub.y] = [[integral of].sub.[gamma]] [n.sub.x] [[sigma].sub.yx][N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.y] [[sigma].sub.yy] [N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.z] [[sigma].sub.yz] [N.sup.i] d[gamma] (13)

[F.sub.z] = [[integral of].sub.[gamma]] [n.sub.x] [[sigma].sub.zx][N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.y] [[sigma].sub.zy] [N.sup.i] d[gamma] + [[integral of].sub.[gamma]] [n.sub.z] [[sigma].sub.zz] [N.sup.i] d[gamma] (14)

In the above equations, [F.sub.i] represents the resultant forces and [N.sup.i] represents the shape functions. Also [N.sub.1], [n.sub.2] and [n.sub.3] are the components of the normal unit vector at the boundary, which for plane [eta]-[zeta] in the normal coordinate system [24] are:

n = [[a.sub.1], [a.sub.2], [a.sub.3]]/[square root of][[a.sup.2].sub.1] + [[a.sup.2].sub.2] + [[a.sup.2].sub.3] (15)

with

[a.sub.1] = [delta]y/[delta][eta] [delta]z/[delta][zeta] - [delta]y/[delta][zeta] [delta]z/[delta][eta] (16)

[a.sub.2] = [delta]x/[delta][zeta] [delta]z/[delta][eta] - [delta]x/[delta][eta] [delta]z/[delta][zeta] (17)

[a.sub.3] = [delta]x/[delta][eta] [delta]y/[delta][zeta] - [delta]x/[delta][zeta] [delta]y/[delta][eta] (18)

A normal stress difference between the entrance and exit surfaces was set as [delta] [[sigma].sub.xx], such that the desired flow rate was obtained. Equations 12, 13 and 14 are applied over the inlet and outlet boundaries. The average pressure rise measured at the apex ports over the length of 42 mm in the 45/5/42 configuration was approximately 650 kPa, and the average pressure rise measured in the side ports was 500 kPa. The reference value of [[sigma].sub.xx] was set zero at the entrance plane and [delta][[sigma].sub.xx] was set to 1035 kPa to match the experimental flow rate of 4.35 Kg/h. the rotation speed was either 60 or 140 RPM.

The problem of the formulation of appropriate boundary conditions at inflow and outflow planes has been a controversial one. Sometimes it is possible to obtain seemingly acceptable solutions for the finite flow field that are, in reality, poor approximations of the real problem [20]. Validation via experimental flow visualization and velocity measurements is the best way to assess model boundary conditions.

2.5 Finite Element Formulation

A three-dimensional finite element method is used to discretize the governing Eqs 1 and 2 in the kneading disc region. Several texts provide the necessary mathematical background of the method applied to fluid mechanics problems [21, 22]. The field variables [velocities and pressures] are approximated by 27-node triquadratic brick elements for velocities and 8 node trilinear brick elements for pressure. For the discretization of the twin screw domain, different densities of 27-node triquadratic brick element meshes were tested. The results shown in this work correspond to a mesh of 2,016 elements and 20,881 nodes, as depicted in Fig. 4.

The finite element code was used to obtain the velocity profiles and the pressure distributions of the flow field according to experimental stagger angles and operating conditions by McCullough and Hilton [1].

3. COMPARISON WITH EXPERIMENTAL PRESSURE PROFILES

Traditional twin screw extrusion process data are limited to pressures and temperatures measured at the die. Measurements are not normally made along the barrel since the best measurement locations are a function of screw configuration. Steep pressure gradients are likely to appear in both the axial and rotation directions. McCullough and Hilton [1] developed an experimental apparatus consisting of a 30 mm diameter TSE using a 2L/D feed adapter, an instrumented 6L/D test barrel section and standard die assembly. The configuration of the apparatus is shown in Fig. 5. The reverse kneading elements at the end of the kneading disc section restricted the flow and raised the pressure within the system. Figure 5a shows an end view of the barrel section from the feed end with the location of the pressure transducer at the top of the apex and the side of the main screw on the right. To obtain the required pressure data, McCullough and Hilton equipped the barrel with five ports evenly spaced every 30 mm along the top and bottom of the center apex region, and three ports along the sides spaced every 60 mm. The TSE assembly was mounted on a movable lathe bed. This allowed the barrel section to be moved back and forth along the screw elements for the various measurements. High-speed data measurements confirmed that the pressure transducers were capable of responding as quickly as 800 MPa/s.

Response time is critical because at 100 rpm (or 600 degrees/s) pressure can fluctuate as much as 8 MPa in 0.05 seconds, or 160 MPa/s. The data collection frequency was set to record pressure data at each degree of screw rotation: 360 records/second at 60 rpm and 800 records per second at 140 rpm.

Results are presented for runs at 60 and 140 rpm using a low density polyethylene. Pressures were measured at the side of the main screw and at the top of the apex, for different arrays of kneading discs, according to Fig. 5b. Of special interest for the comparisons with simulation data are two of the configurations depicted in Fig. 5b. The configuration of forward kneading discs are 45/5/20 (where 45 is the stagger angle, 5 is the number of discs and 20 stands for a total disc section length of 20 mm) and 45/5/42. The pressure measurements are plotted as a function of the amount of rotation. The axial pressure drop plots show average, maximum and minimum values over one cycle of rotation for apex and side measurements. For the simulations, the average side pressure drop was used as an initial guess for the normal stress difference imposed as boundary condition and then adjusted to obtain the correct flow rate. The average side pressure was used since it is the most representative of the average behavior in the system. The pressure fluctuations in the apex region are more intense than at the side because of the interaction of the two screws at the intermeshing region. The flow rate was kept at 4.35 kg/hr to match experimental conditions The complete set of operating conditions and material properties (1) are presented in Tables 2 and 3 and Fig. 6.

A total of 7 steady state solutions with increments of 12.86[degrees] are used to complete one quarter of a revolution. Because of the symmetry of the system, calculation of only one quarter of a revolution is necessary to describe the complete cycle. The complete set of solutions is obtained by means of spatial transformation to rotate the first quarter geometry 180[degrees] around the x-axis to obtain the second quarter. Once the solution for a half revolution is obtained, symmetry applies again making the second half revolution identical to the first half. An isometric view of the combination of kneading blocks used for the experimental work is shown in Fig. 7.

The numerical determination of pressures at the side and apex is done by interpolating the values at the locations equivalent to the location of the pressure transducers in the experimental setup. The reported experimental pressure value is centered at a point, measured over a small area. Nine calculated pressure values are interpolated in a cross pattern as shown in Fig. 8 to cover a 1-mm-diameter circle, which corresponds to the area exposed to the pressure transducer in the experimental setup. An average of the nine points is then taken to account for the pressure on the 1-mm-diameter circle. The pressure values are interpolated from the nodal values previously obtained through the FEM calculations. Once the coordinate of the point where the pressure needs to be calculated is determined, an algorithm searches for the element of the mesh that contains that coordinate. The pressure is then interpolated from the nodal values using the element's interpolation functions. For a full revolution, the numerical pr ocedure of determining the pressure at the apex and side ports is repeated 28 times for every disc.

For the simulation of the 45/5/42 configuration the geometry depicted in Fig. 5 representing 5 pairs of discs was used. For the 45/5/20 configuration a system with two pairs of discs was used. Two disks were used for the 45/5/20 configuration to save computer time, given the limited computational resources. However, from Fig. 5b, there is an important effect of disk 4 of the 45/5/20 block on the circumferential pressure profile of disk 5, due to their proximity. Disk 1 of the 45/5/42 is aligned with disk 5 of the 45/5/20 block, and disk 2 of the 45/5/42 block does not have a major effect. Experimental and simulation results have a better match for this particular case, because of the distance of disk 2 of 45/5/42 from disk 5 of 45/5/20. However, a more complete analysis would include more disks in the 45/5/20 configuration.

4. DISCUSSION

Figures 9 and 10 show the comparison between simulation and experimental variation in pressure for a complete cycle of disc 5 in the 45/5/20 configuration. At this location, the pressure fluctuations are the largest in the system. The experimental peaks for the apex port reach 4250 kPa whereas the minimum pressure is 1750 kPa. Figure 10 shows the pressure profiles at the side port, where the amplitude of the fluctuations is less than for the apex port. Figures 11 to 13 present experimental and simulation results for the individual discs of the 45/5/42 configuration at the apex port. Figures 14 to 16 show the experimental and simulation results at the side port of the 45/5/42 configuration.

An increasing pressure in the axial direction has been reported [1] for the configuration of forward kneading blocks, as predicted by the simulations. Although the staggering of the discs reduces the conveying capacity of kneading blocks, they still maintain some ability to move material forward. The forward conveying is produced by the dragging effect whereas the axial pressure induces a back flow component that ensures a fully filled channel and hence better mixing performance. The normal stress boundary condition determines whether there is a positive or negative conveying effect for specific pressure drop and rotation speed conditions.

The motion of the rotors causes the fluctuating pressures observed experimentally and predicted with the model. Different mechanisms cause the pressure fluctuations in the side and apex. At the side observation point, the advancing tip of the discs creates a converging region between the disc and the barrel causing high-pressure peaks. The relative motion of the trailing side of the disc to the barrel wall produces the low pressures, as illustrated in Fig. 17. The high pressures at the apex originate from the encounter of the tips of the two discs. A small chamber is created and then reduced in size until it completely disappears. squeezing the fluid and creating the high peaks of pressure. Immediately after the chamber disappears, a small chamber appears again and increases in size until it communicates with the wide channels at both sides of the chamber, creating the low pressure values, as illustrated in Fig. 18.

Pressure valleys at one disc coincide with the highest pressure peaks in the subsequent disc, generating the experimentally observed back flows at the rear side of the discs. This effect is very much affected by the stagger angle as well as the thickness of the discs. It is expected that thicker discs would generate wider valleys of pressure and therefore greater back flows. The 45[degrees] stagger angle kneading blocks have the particular property of the coincidence of the highest pressure peaks on the advancing side of the disc's tip with the lowest pressure valleys on the rear side of the contiguous disc's tip. Based on this observation and from the model calculations, the phenomena observed by McCullough and Hilton [1] can be explained. Figures 9 and 10 show the pressure behavior as a function of the angle of rotation for disc 5 in the 45/5/20 configuration. Disc 5 is adjacent and aligned with the first disc of the 45/5/42 configuration (Fig. 5b). The resulting effect is that of a long disc by adding the 4 mm thickness of the 45/5/20 disc 5 to the 8.2 mm of the 45/5/42 disc 1, for 12.2 mm total thickness. Simulations reveal that the pressure peaks for a single pair of discs, under similar flow rate and rotational speed to those used in the experimental work, are higher than the peaks produced in discs preceded and followed by another disc as occurs in real kneading blocks. This behavior has its explanation in the fact that the peaks and valleys of pressure for one disc are dampened by the effect of discs behind and ahead of it. In the case of disc 5 of the 45/5/20, the subsequent downstream staggered disc is at a distance of 8.22 mm, thus reducing its effect over the pressure profile of the preceding disc. This explains the experimental observation of an increase in the pressure fluctuations from the peaks to the valleys when approaching disc 5. It also explains why the amount of fluctuation decreases rapidly once the 45/5/42 configuration is reached. The intensity of pressure fluctuations is affected by the reverse elements located at the end of the 45/5/42 configuration that actually holds up the polymer melt and creates an almost uniform pressure distribution over the cross section of disc 5 of the 45/5/42 block. The simulation of the 45/5/42 configuration (Fig. 5b) used 5 pairs of discs. For the 45/5/20 configuration, two pairs of discs were used in the simulation. Disc 5 of the 45/5/20 block is aligned with the first disc of the 45/5/42 block, thus its local region is different from the rest of the discs in the same configuration. Comparisons of the model with experimental data on the 45/5/20 configuration (Figs. 9 and 10) show very good correlation, and together with the simulations for the 45/5/42 configuration help to explain the pressure behavior observed experimentally. Simulations for the 45/5/42 show also good correlation for all discs but disc 5, which is expected since it is adjacent to the reverse kneading blocks.

Although the model predictions are reasonably close to experimental pressure data, there is a consistent difference in the pressure peaks. When analyzing the model, there are two main assumptions made among those discussed in Sections 2 and 3: quasi-steady state and isothermal behavior. The quasi steady state assumption is said to apply for systems where Re (Reynold number, Re = [rho]VD/[micro]) [less than][less than]1 and Sr (Strouhal number Sr = [omega]D/V) [less than][less than]1. The assumption of Re [less than][less than]1 for polymeric flows is valid due to the very high viscosity of polymer melts and hence acceleration terms are neglected. The assumption of Sr [less than][less than]1 implies that the time scale of the perturbation, i.e. a periodic change in boundary conditions such as a pulsating flow or a moving boundary, must be much larger than the time scale of the diffusion of momentum. i.e. t[greater than][greater than][D.sup.2]/[nu]. This assumption is valid for polymeric flows that demonstrate ve ry low values of [D.sup.2]/[nu]. There are very large pressure gradients generated in a very small time interval due to the sudden transition from the narrow gap of the disc-barrel clearance to the wider channel that follows. These pressure peaks represent extremely small zones within the flow field that may not affect flow quantities for purposes such as determining mixing performance from numerically calculated flow fields. Inclusion of transient terms and an algorithm to integrate the equations in time may capture the sharp pressure transitions accurately.

The system is assumed to be isothermal. In real systems, it is expected that zones of high shear would experience higher temperatures due to viscous dissipation and therefore lower values of viscosity than estimated from the isothermal model. It is well known that the zones where the pressure peaks occur are also the zones with the highest shear rates in the system. The isothermal flow assumption could be a reason for overestimation of the pressure peaks in the system.

An additional reason for the discrepancies between experimental and simulation pressures could be the model for the dependence of viscosity on the magnitude of the rate-of-strain tensor [gamma] (see Eqs 4 to 7). The highest peaks of pressure occur at the narrow gaps where the [gamma] values are the highest. The Carreau model fits the data for a range of [gamma] between 0 and 5000 [s.sup.-1] according to Fig. 6. Calculations have demonstrated that at the narrow gaps the values of [gamma] can be greater that 5000 [s.sup.-1]. Over-prediction of the viscosity at high values of [gamma] could produce an overestimation of the pressure profiles. An alternative to the present model, would be a piece-wise fit of the Carreau model to specified ranges of shear rate to capture the apparently lower power law index at shear rates in the 1000-10,000 [s.sup.-1] range.

The results are also in qualitative agreement with circumferential pressure measurements presented by Kiani and Heidemeyer [23] and calculations made by Szydlowski et al. [4]. Physical limitations in the pressure transducers--i.e., speed of response, response to a non-uniform pressure in the face of the transducer, noise, etc.--may also introduce deviations from the prediction made by the model. Pressure transducers have limitations in the measurement of steep gradients over small areas since they are usually designed to measure constant pressures over the area of the probe. Limitations in the speed of response of the transducer could produce a "dampening" of the signal at those locations where the gradients are the steepest.

5. CONCLUSIONS

A quasi-steady state 3-D finite element model has been implemented and compared against experimental data provided by McCullough and Hilton [1]. For the purpose of comparison, simulation results were analyzed over small areas equivalent to that of the pressure transducers, since numerical interpolation determines the pressure value at a point. Differences are noticed at the pressure peaks produced by the tips of the kneading discs "sweeping" the surface of the barrel, whereas for the other regions the comparison is reasonably good. Sources of error to explain the differences between simulation and experiment include: mesh density, isothermal assumption, fit of the constitutive equation at high shear rates, and pressure transducer response time.

ACKNOWLEDGMENTS

The authors wish to acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada and Investigacion y Desarrollo C.A. (Venezuela).

(1.) V. L. Bravo on leave from Investigacion y Desarrollo C.A. (INDESCA), P.O. Box 10319. Complejo Petroquimico El Tablazo. Maracalbo, 4001. Venezuela.

(*.) Author to whom correspondence should be addressed

(email: hrymak@mcmaster.ca, fax: (905) 521-1350).

6. REFERENCES

(1.) T. W. McCullough and B. T. Hilton, SPE ANTEC Technical Papers, 3372 (1993).

(2.) A. Kaplan and Z. Tadmor, Polym. Eng. Sci., 14, 58 (1974).

(3.) Z. Tadmor, E. Broyer, and C. Gutfinger, Polym. Eng. Sci., 14, 660 (1974].

(4.) W. Szydlowski and J. L. White, Intern. Polym. Processing, 2, 142 (1987).

(5.) M. L. Booy, Polym. Eng. Sci., 18, 3 (1978).

(6.) M. L. Booy, Polym. Eng. Sci., 20, 1220 (1980).

(7.) A. D. Gotsis, Z. Ji, and D. M. Kaylon, SPE ANTEC Technical Papers. 139 (1990).

(8.) A. Lawal and D. M. Kalyon, Polym. Eng. Sci., 35, 1325 (1995).

(9.) Y. Wang. J. L. White, and W. Szydlowski, Intern. Polym. Proc., 4, 262 (1989)

(10.) A. Kiani, E. W. Grald, and S. Subbiah, "Spectral Element Analysis of the Mixing Characteristics of Twin Screw Extruders," Session on Numerical Simulation of Mixing Phenomena. AIChE Annual Meeting (1992).

(11.) A. Kiani and H. J. Samann, SPE ANTEC Technical Papers, 2578 (1993).

(12.) A. T. Patera, J. Computational Physics, 54, 468 (1984).

(13.) H. H. Yang, Flow Field Analysis of Batch and Continuous Mixing Equipment, PhD thesis. Case-Western Reserve University, Cleveland, Ohio (1993).

(14.) I. Manas-Zloczower and J. J. Cheng, Polym. Eng. Sci., 29, 1059 (1989).

(15.) T. Kajiwara, Y. Nagashima, Y. Nakano, and K. Funatsu, Polym. Eng. Sci., 36, 2142 (1996)

(16.) C. Rauwendaal, Polym. Eng. Sci., 21, 1092 (1981).

(17.) J. Christiano and M. Lindenfelzer, SPE ANTEC Technical Papers, 78 (1997).

(18.) R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena. John Wiley (1984).

(19.) V. Bravo, 3D Simulation of Quasi-Steady Flow Field in Twin Screw Extruders, M. Eng. thesis, McMaster University. Hamilton, Ontario, Canada (1995).

(20.) J. S. Vrentas, C. M. Vrentas, and H. C. Ling, J Non-Newt. Fluid Mech., 19, 1 (1985).

(21.) G. Dhatt and G. Touzot, The Finite Element Method Displayed, John Wiley and Sons, New York (1984).

(22.) J. N. Reddy, Introduction to the Finite Element Method, McGraw-Hill, New York (1984).

(23.) A. Kiani and P. Heidemeyer, SPE ANTEC Technical Papers, 94 (1997).

Technical Data, ZSK 30 Extruder. Barrel bore diameter ([D.sub.b]) 30.85 mm Screw outside diameter ([D.sub.o]) 30.7 mm Screw root diameter ([D.sub.l]) 21.3 mm Flight depth (t) 4.7 mm Center distance of screws ([C.sub.L]) 26.2 mm Clearance between discs (G) 0.2 mm Clearance disc-barrel (B) 0.075 mm Operating Conditions. Flow rate 4.35 Kg/hr Rpm 60/140 Kneading blocks 45/5/20, 45/5/42 Average melt temperature 223[degrees]C Location of ports Apex and side Material Properties and Rheological Parameters. Material Properties Material LDPE 6411 (Dow Chemical) Melt Density 740 Kg/[m.sup.3] Rheological parameters for Carreau model [[micro].sub.o] 1,290 Pa.s [lambda] 0.112 s N 0.559

Printer friendly Cite/link Email Feedback | |

Author: | BRAVO, V.L.; HRYMAK, A.N.; WRIGHT, J.D. |
---|---|

Publication: | Polymer Engineering and Science |

Article Type: | Statistical Data Included |

Geographic Code: | 1USA |

Date: | Feb 1, 2000 |

Words: | 6736 |

Previous Article: | Dynamic Slip and Nonlinear Viscoelasticity. |

Next Article: | Preparation and Properties of Thermoplastic Starch-Polyester Laminate Sheets by Coextrusion. |

Topics: |