3-D Numerical Simulations of Nonisothermal Flow in Co-Rotating Twin Screw Extruders.
KAZUMORI FUNATSU [*]
Three-dimensional nonisothermal flow simulations in the kneading disc regions of co-rotating twin screw extruders were performed using a finite element method. The standard Galerkin method and penalty function scheme were applied to the flow field. The streamline-upwind/Petrov-Galerkin scheme was used in the temperature field to reduce numerical oscillation. The simulations were carried out under the operational conditions of The Japan Steel Works TEX3O machine for various rotational speeds. The configuration was ten 2-lobe kneading discs with a 90[degrees] stagger angle. Experimental observations were also performed to validate the numerical simulations under the same operational conditions. The pressure in front of the tip in the rotation direction was higher than behind the tip, and the region behind the tip sometimes had a negative value. Since variation of the pressure gradient in the axial direction causes forward and backward flows in the disc gap regions, the disc gap regions play an important role f or mixing. The temperature becomes higher with increasing rotation speed due to high viscous dissipation. A high temperature was observed on the disc surface, in the disc gap, and in the intermeshing regions. The numerical results of pressure profiles with the rotation and the temperature in the axial direction were in good agreement with the experimental observations.
Twin screw extruders are widely used in the plastics industry for distributing and dispersing of polymers and additives. Many kinds of twin screw extruders have been presented, and they are mainly classified by the rotation direction of each screw and the degree of intermeshing. In these machines, intermeshing co-rotating twin screw extruders are the most popular in the polymer industry and consist of two types of screw elements. One is a full-flight screw, and it is used for conveying polymer melt from the feed section towards the die exit. Another type of screw element is a kneading disc, and it plays an important role in the dispersive and distributive mixing of polymeric materials.
Since it is desirable to characterize the polymer flow in twin screw extruders theoretically, many numerical studies have been developed. However, sufficient results have not been obtained due to the complicated 3-D geometry of the kneading discs. One numerical approach uses the approximation technique for a low computational cost. Szydlowski et al.  simulated the flow in the kneading disc regions of co-rotating twin screw extruders by means of a Flow Analysis Network (FAN). The FAN technique uses lubrication approximation and neglects the flow in the radial direction. This 2-D assumption caused several geometrical problems. The change of intermeshing geometry with screw rotation was improved by Szydlowski et al.  and the curvature effect was considered by Sebastian et al. . Another numerical approach uses the Finite Element Method (FEM). Yang et al.  and Li et al.  developed the simulation technique of isothermal flow of mixing element regions using a commercial FEM code. Kajiwara et al. [6, 7] also developed a FEM program of isothermal flow in co-rotating twin screw extruders. Since the FEM technique has no limitations of geometry, it was considered to be suitable for the kneading disc regions treated in this study.
In this study, we developed a nonisothermal simulation program in the kneading disc regions of twin screw extruders using the FEM technique. The simulations were carried out under the operational conditions of a Japan Steel Works TEX30 machine at various rotational speeds. The numerical results, such as pressure, velocity and temperature distribution, were precisely examined. Experimental observations of pressure profiles with rotation and temperature in the axial direction were also carried out and compared with numerical results.
In the case of twin screw extruders, the flow domain changes with screw rotation. This moving boundary problem becomes the major difficulty in numerical analysis. To solve this, Yang et al.  performed a number of quasi-steady state analyses of the sequential geometries to represent dynamic motion. Since the Reynolds number of polymer flow is very small, the isothermal flow can be approximated as quasisteady state. Since the thermal conductivity of the polymer melt is small in the case of nonisothermal flow, the temperature distribution Is dominated by viscous heating and convection. The Peclet number and the Brinkman number are on the order of [10.sup.2] at a rotational speed of 200 rpm. If the temperature distribution is mainly dominated by convection, the quasi-steady state analysis would not be acceptable. However, in this study we assumed that the temperature distribution was mainly dominated by transient viscous heating. Hence, the quasi-steady state analysis could be applied. Although this assumptio n is not appropriate at a high rotational speed, the obtained results can be useful for the relative evaluation of the flow field according to the comparison with the following experiments. We also verified by experimental observation that the flow domain was fully filled with polymer melt.
The model is based on pure viscous non-Newtonian and incompressible flow with the stress tensor [tau], pressure p, and velocity vector v. The continuity and momentum equations are
[nable] * v = 0 (1)
- [nabla]p + [nabla]*tau] = 0. (2)
The energy equation with density p, temperature [rho], temperature T, heat capacity [C.sub.p], and thermal conductivity k is
[rho][C.sub.p]v* [nabla]T = k[[nabla].sup.2]T + [tau] : [nabla]v. (3)
The constitutive equation and its temperature dependency are represented by the Carreau model and Arrhenius' law.
[tau] = 2[eta]D (4)
[eta] = H(T)F([gamma] H(T)), [gamma] = [square root of]2[[pi].sub.D] (5)
F = [[eta].sub.0][[1 + [[[lambda].sup.2].sub.c](2[[pi].sub.D])].sup.(n-1)/2] (6)
H(T) = exp[-[beta](T - [T.sub.a])] (7)
where D is the rate deformation tensor, [[eta].sub.0] the zero shear rate viscosity, [[pi].sub.D] the second invariant of the rate deformation tensor, [[lambda].sub.c] and n the Carreau model's parameters, [beta] the Arrhenius law's parameter, and T[alpha] the reference temperature, respectively.
The standard FEM Galerkin method and penalty function scheme were used in deriving the flow field. The SUPG [8, 9] scheme was applied to the temperature field to obtain oscillation-free solutions.
We performed the analysis in a flow field and a temperature field separately to reduce the computer memory and get convergence via iterative scheme. The flow diagram of this simulation is shown in Fig. 1. The matrix solution methods were the skyline method for the flow field and the conjugate gradient method for the temperature field, respectively.
The finite element design of flow domain in The Japan Steel Works TEX30 disc and barrel surfaces are shown in Fig. 2. The configuration consists of ten 2-lobe kneading discs with a 90[degrees] stagger. The disc dimensions and the process conditions are listed in Tables 1 and 2 in detail. The boundary conditions of the flow and the heat analysis are listed in Tables 3 and 4, respectively.
In order to represent the dynamic motion, we choose six different geometries in a sequential manner with 15[degrees] increments, due to the symmetry of the discs. The geometries representing a quarter of the rotation cycle are shown in Fig. 3 and are labeled by the angle [alpha] between the X-axis and the left disc tip.
The material used in this study was polypropylene produced by Mitsubishi Chemical Corp., and its viscosity was measured on a capillary rheometer Goettfert Rheograph 2002 at a high shear rate and on a mechanical spectrometer Rheometric RMS-800 at a low shear rate. The curve fitting results by the Carreau model are shown in Fig. 4 and its parameters are listed in Table 5. The other parameters in Table 5 are data typical of polypropylene.
In the numerical simulations, the number of nodes was 22,561 and the degree of freedom was 60,843 in the flow field and 20,281 in temperature field. The required computer memory and time were about 480 MB and 90 minutes by a workstation HP-J2240 without parallel computing.
The measurements of pressure and temperature were also carried out under the same configuration and the same operational conditions as the numerical simulations. Figure 5 shows the barrel section equipped with several ports for pressure and temperature measurements. Four ports are evenly spaced every 18.5 mm for temperature measurements in the nip (intermeshing) region. A port for pressure measurements is located at the same axial direction as the first temperature port in the translation region (away from the nip). The pressure gauge, K. K. Dynisco NP462-1/22CK-15/45, and the thermocouples, Okazaki Manufacturing Co. Okazaki Aeropak Sheathed Thermocouples T90 Type K, were placed in the port and online measurements were performed.
The position of the cross section used for contour plots of the numerical results is also shown in Fig. 5.
RESULTS AND DISCUSSION
The pressure distributions at the rotational speed 200 rpm on different cross sections (see Fig. 5) are shown in Fig. 6. Figure 6A shows the upper stream disc gap cross section A-A', Fig. 6B shows the disc center cross section B-B' and Fig. 6C shows the downstream disc gap cross section C-C'. The pressure in front of the tip in the rotational direction is shown to be higher than that behind the tip in all cases, and the region behind the right upper tip in Fig. 6B is shown to be negative. The pressure gradient in the Z direction changes with its Z position. For example in the region behind the upper right tip in the rotational direction, the pressure gradient (A)-(B) shows a negative value, and (B)-(C) shows positive. On the other hand, in the region in front of the lower right tip in the rotational direction, the pressure gradient (A)-(B) shows a positive value and (B)-(C) shows negative. These pressure gradients are caused by forward and backward displacement of polymer melts in the Z direction due to disc rotation. The comparison between numerical and experimental results of pressure profiles with disc rotation at the rotational speed 200 rpm is shown in Fig. 7. It shows good agreement in terms of the pressure value and oscillation period. The point of pressure measurement is shown in Fig. 5.
Figure 8 shows the velocity vectors at the rotational speed 200 rpm on a B-B' cross section, and the maximum velocity is 0.37m/sec. Figure 9 shows the Z component of velocity [V.sub.z] on a different cross section in the same position as in Figs. 5 and 6. The velocity component [V.sub.z] in (B) is very small in comparison with (A) and (C). A large amount of forward and backward flow in the Z direction was observed in (A) and (C). For example, the value of [V.sub.z] is positive at the region behind the upper right tip in the rotational direction and negative at the region in front of the lower right tip in (A). On the other hand, in (C) the reverse distribution is observed. It is caused by the variation of the pressure gradient in the Z direction as mentioned before. We understand the disc gap regions (A) and (C) play an important role for mixing due to the large amount of forward and backward flow.
The temperature distributions on the B-B cross section at a different rotational degree [alpha] at a rotational speed of 200 rpm are shown in Fig. 10. The barrel surface temperature is lower than the disc surface temperature, and the temperature in the intermeshing region is highest in all cases. This is caused by the cooling effects from the temperature-controlled barrel and the adiabatic boundary condition on the disc surface. Figure 11 shows the temperature distributions on disc and barrel surfaces at different rotational speeds. The upper left region of each domain is removed for visualization. The temperature becomes higher with increasing rotational speed, and the intermeshing and the disc gap regions show higher values in all cases. The comparison between numerical and experimental results of temperature in the Z direction at different disc rotational speeds is shown in Fig. 12. In this figure, the temperature of numerical results is the average value of each rotational degree at the position shown in Fig. 5. The numerical results show good agreement with the experimental results at low rotational speeds of 100 and 200 rpm. At the high rotational speed of 400 rpm, however, the numerical result somewhat overestimates the experimental result. It is considered to be the effect of the quasi-steady state assumption, which causes small thermal convection. Another reason is overestimation of viscous dissipation due to a wall slip .
SUMMARY AND CONCLUSIONS
In this paper, a 3-D nonisothermal flow program using FEM was described. The simulations were performed in the kneading disc regions of twin screw extruders at various rotational speeds. Experimental observations of pressure and temperature were also carried out. The results were precisely examined and the following conclusions were obtained:
(1) The pressure in front of the tip in the rotational direction was shown to be higher than that behind the tip, and the region behind the tip sometimes was shown to have a negative value. The pressure gradient in the Z direction changed with its Z position because of disc rotation.
(2) The forward and backward flows were observed in the disc gap regions due to the variation of the pressure gradient. The disc gap regions are considered to play an important role for mixing.
(3) The temperature becomes higher with increasing rotational speed and the values on the disc surface, in the disc gap, and in the intermeshing regions are higher in all cases.
(4) The numerical results show good agreement with experimental observations in terms of the pressure profile with disc rotation and the temperature distribution in the Z direction.
We extend our thanks to Mr. Teruo Amaiwa and Mr. Tsukasa Sato of Mitsubishi Chemical Corp.'s Yokkaichi plant for useful discussions and experimental support.
(1.) 0n leave from Mitsubishi Chemical Corp., Yokkaichi plant 1. Toho-cho, Yokkaichi 510-8530. Japan
(*.) To whom correspondence should he addressed.
(1.) W. Szydlowski, R. Brzoskowski, and J. L. White, Intern. Polym. Proc., 1. 207 (1987).
(2.) W. Szydlowski and J. L. White, Intern. Polym. Proc., 2, 142 (1988].
(3.) D. H. Sebastian and R. Rakos, SPE ANTEC Tech. Papers, 36, 135 (1990).
(4.) H-H. Yang and I. Manas-Zloczower, Polym. Eng. Sci., 32, 1411 (1992).
(5.) T. Li and I. Manas-Zloczower, sPolym. Proc., 10, 314 (1995).
(6.) T. Kajiwara, Y. Nakano, S. Ninomiya, and K. Funatsu, Seikei-Kakou, 5, 557 (1993), (in Japanese].
(7.) T. Kajiwara, Y. Nagashima, Y. Nakano, and K. Funatsu, Polym. Eng. Sci., 36, 2142 (1996).
(8.) A. N. Brooks and T. J. R. Hughes, Camp. Meth. Appl. Mech. Eng., 32, 199 (1982).
(9.) J. M. Marchal and M. J. Crochet, J. Non-Newt. Fluid Mech., 26, 77 (1987).
(10.) X. Yang, H. Ishida, and S. Q. Wang, J. Rhoel., 42, 63 (1998).
Dimensions of 2-Lobe Kneading Disc Regions. Barrel diameter D [mm] 30.0 Screw tip diameter [D.sub.S] [mm] 29.2 Screw root diameter [D.sub.R] [mm] 21.0 Centerline distance [C.sub.L] [mm] 26.0 Disk width W [mm] 7.8 Disk gap [[delta].sub.s] [mm] 1.5 Stagger angle [psi] [deg] 90 Number of Discs [N.sub.D] 10 Process Conditions. Barrel temperature [T.sub.B][[degrees]C] 200 Rotational speed N [rpm] 100, 200, 400 Flow rate Q [kg/h] 15 Flow Boundary Conditions. Inlet cross section Constant flow rate (15 kg/h) Barrel inner surface No slip Disc surface Tangential velocity by screw rotation Outlet cross section Outflow Temperature Boundary Conditions. Inlet cross section Constant temperature 202, 208, 216[degrees]C at 100, 200, 400 rpm Barrel inner surface Heat flux: q = h(T-Ta) h = 3000 [W/[m.sup.2]K] Ta = Inlet cross section temperature Disc surface Adiabatic Outlet cross section Adiabatic (dT/dz = 0) Material Data of Polypropylene. Carreau model param. [lambda] 2.15 Carreau model param. n 0.38 Zero shear rate viscosity [[eta].sub.0]] [Pa.s] 26470.0 Arrhenius' law param. [beta] [1/K] 0.02 Reference temp. [T.sub.[alpha]] [K] 473.0 Density [rho] [kg/[m.sup.3]] 735.0 Specific heat [C.sub.p] [J/kgK] 2100.0 Thermal conductivity k [W/mK] 0.15
|Printer friendly Cite/link Email Feedback|
|Author:||ISHIKAWA, TAKESHI; KIHARA, SHIN-ICHI; FUNATSU, KAZUMORI|
|Publication:||Polymer Engineering and Science|
|Article Type:||Brief Article|
|Date:||Feb 1, 2000|
|Previous Article:||Modeling of Peroxide Initiated Controlled Degradation of Polypropylene in a Twin Screw Extruder.|
|Next Article:||Numerical Simulation and Experimental Verification of Nonisothermal Flow in Counter-Rotating Nonintermeshing Continuous Mixers.|