Numerical Modeling of Thermal Flows in Entrance Channels for Polymer Extrusion: A Parametric Study

: Flow distribution channels in extrusion dies are typically designed to assure uniform ﬂuid velocity, pressure and temperature in the outlets. To ensure this uniformity, it is desirable to have the ﬂuid melt to reach a steady state temperature in the entrance channel before entering the die body. This paper numerically investigates the temperature distribution of the ﬂuid melt in the entrance channel. Analytical solutions of the velocity and ﬁnite element solutions of temperature distribution in Poiseuille ﬂows of polypropylene melt with the Casson rheology model were derived and presented. In the velocity solution, the critical point that separates the core and the remaining parts in the ﬂow was calculated by using the inlet ﬂow rate and the yield stress in the Casson model. The velocity distribution was then substituted into the convective heat equation for temperature distribution simulations. A ﬁnite difference scheme was used to obtain the temperature distribution proﬁles along the ﬂow direction in a parallel-plate, while the ﬁnite element model was used to model the ﬂow temperature in circular tubes. The main outcome is the parametric analyses of the effect of various parameters such as radius, wall temperature, inlet temperature, and pressure drop to the optimal length of the channels required for the ﬂow temperature to reach the steady state.


Introduction
It is well known that polypropylene can be used to produce low-cost, lightweight, transparent, and flexible thin sheets and films for a variety of applications in industries [1]. The key component in the production of these sheets and films is the extrusion die designs [2,3]. Flow distribution channels in the extrusion dies should be designed in such a way to assure uniform fluid velocity, pressure, temperature in the outlets with uniform thickness [4,5]. If the fluid flow is rapid at the centerline or if the uneven temperature distribution occurs then these scenarios result in extrusion product deformations and surface cracks [6]. Therefore, understanding the Poiseuille flows of polypropylene melt is fundamental for the design of extrusion dies for production of high quality PP sheets and PP films [7,8].
Poiseuille flows of non-Newtonian fluids are well studied in many rheology models [9,10]. For most of the rheology models, analytical and numerical velocity profiles have been available in standard textbooks and literature [11]. The Casson model has often been used to model Non-Newtonian fluids in food industries and it has been used for modeling blood flows [12,13]. Moreover, it can also be used for modeling the polymer flows, particularly, for polypropylene melts.
Although the velocity profile is well-known in literature, the connection of the profile with the inlet flow rate and the core size is not explicitly explained [14,15]. When the polymer melt enters an entrance pipe of the die, the melt is at the melting temperature but the die body is kept at a higher temperature anticipating that the melt temperature will rise due to Arrhenius law and the shear rate of the melt at the wall [16,17]. Similar research work has been done by Wei and Luo, where they numerically analyze a power law heat transfer problem of polymer melt flows in a tube with constant ambient temperature [18]. However, in this paper, we apply the Casson fluid model to solve the temperature distribution of the PP melt flow. In general, it is commonly used to describe the blood, chocolate and ketchup flows [19][20][21]. However, some studies mention that the Casson model is applied in polymer industries by considering it for the rheology analysis of the polymer melt flows [22,23]. Moreover, the studies of Lungu et al. [24] suggest that there is a good fit of the Casson model for the polymer melt flows' experimental data. In addition, Matveenko and Kirsanov [25] showed that the generalized Casson rheological model describes the behavior of polymer melt flows better than the well-known and commonly used power-law model.
In general, polymer extrusion die has different types of design used in the industry, for instance, coat hanger or T-slot dies [26]. Since our focus is the extrusion die, the parts of the die we are analyzing are the circular tube that often appears to be the inlet part and the parallel-plate structure that is sometimes used in the manifold. In this paper, the finite difference was used to solve the convective heat equation because of the simplicity of the model. However, for the flow temperature simulation in the circular tube, the finite difference scheme was no longer applicable due to the singularity occurring at r = 0 in the solution of the nonlinear partial differential equation. Therefore, a finite element model was found to be appropriate since it additionally applies the derivative boundary condition.
Eventually, it is desirable to have the fluid melt to reach a steady state temperature before entering into the die body for further flow distribution. It is the purpose of this work to estimate the distance from the entrance along the pipe at which the steady state is reached. The entrance pipe is usually circular, so the axisymmetric nonlinear finite element model was developed for the simulation of the temperature when the fluid enters the die and gets distributed into die pre-land, which is in the form of parallel-plates ( Figure 1). For further steady state temperature distribution inside the slits, the finite difference scheme was also developed for the simulation. The results can be used as the first steps leading to calculations of the corresponding fluid velocity and temperature distributions of the downstream extrusion die channel to the manifold and pre-land of a flat die for PP melt flows.

Mathematical Model of Temperature Distribution of the Casson Fluid Flow
Considering PP flows between two parallel plates shown Figure 1 below, we will derive the velocity profile by determining the core velocity u c , the separation point y c between the core and the non-core flows, and the velocity of the non-core portion u, and finally the temperature distribution T along the two-dimensional domain in the following subsections.

Derivation of the Velocity Profile
For Poiseuille flows, the Casson model has the following form where du dy is the shear rate, τ is the shear stress, τ 0 is the Casson yield stress; and η c is a Casson viscosity. The velocity for the case |τ| ≤ τ 0 is constant. Now, we solve for the case |τ| ≥ τ 0 . From the Casson model and balance of forces, the velocity u can be derived as follows where -dP dz is a uniform pressure gradient, L . The integration of both sides from the core boundary y c to the wall y max , we have The integration of shear rate will give the velocity in the interval y c ≤ y ≤ y max .
At the critical point y c , the velocity core can be obtained (constant velocity about center line) In order to find the critical point where du dy = 0 and τ = τ 0 Q = y max y c 2udy + 2u core y c We solve the previous equation, and the equation below is numerically solved to find the critical y c We numerically calculate the separation point y c from Equation (8), which is 0.18 mm by using the constants from [27][28][29], i.e., η c = 237.2 Pa.s, τ 0 = 11786.1 Pa, y max = 0.0018 m, g = 12 MPa. The plot of the velocity profile in Figure 2 is obtained as an example.

Temperature Distribution of the Flow in Parallel-Plate Channel
Convective heat equation for the fluid is given by the balance of energy equation with the Dirichlet boundary condition T inlet = 403.5 K and T wall = 433.5 K where ρ is a fluid density, C p is a specific heat capacity, and k is a fluid thermal conductivity. For the simulation, Equations (3)-(5) are used in solving the heat Equation (9) with constant values shown in [30]. We will use the Crank-Nickolson method with the central difference formula with where i = 0 and i = m + 1 are for boundary points, and i = 1, 2, . . . , m are for interior points, and l = 0, 1, 2, . . . , m represent the steps along z [31]. Next, we approximate the first and second order derivatives with the centered formula Substituting the above approximations into the heat Equation (9) yields to Following the regular finite difference procedure, we represent the above expression as the system of algebraic equations for interior points i = 1, 2, . . . , n and we obtain the equation in matrix form A{T}={R} where A is n × n matrix and {T}, {R} are n × 1 vectors as defined by (12), (13), and T 0 , T n+1 are temperature values at boundary points.
With z as a direction of the flow and y as a vertical position of the fluid, the temperature distribution of the flow in parallel-plate can be calculated by using the Thomas algorithm, which is a highly efficient method for solving the matrix equations with tridiagonal form [31].

Temperature Distribution of the Flow in Circular Tube Channel
The convective heat equation in polar coordinates is given by the balance of the energy equation where the constants are the same as in Equation (9). The boundary conditions are T(r, 0) = T inlet , T(r 0 , z) = T wall , ∂T(0,z) ∂r = 0 and ∂T(r,z) ∂z = 0 as z − → ∞. According to Wei and Luo [18], T(r 0 , z) = T wall will be approximated by the mixed boundary condition −k ∂T(r,z) ∂r = h(T(r, z) − T wall ). The dimensionless parameters are as follows where z 0 is a maximum distance along z direction and r 0 is a maximum radius of the tube. The non-dimensional form of the Equation (14) is where D = . For the circular tube flow, the dimensional velocity with dimensionless parameters is and the corresponding equation of change in velocity is The equation of yield stress defined by the Casson model is

Finite Element Model
We approximate T byT = W T T where T = [T i T j T k ] T and W = [W i W j W k ] T is the vector of three linear shape functions of (R, Z) [32]. Following the Galerkin approach, we have the weak integral form where D e is the element domain. The second volume integral in (19) can be reduced to a surface integral by the Divergence Theorem and the derivative boundary is given by the normal flux, so we have where G e is the boundary domain. Before we construct the local stiffness matrices and load vectors for (19), let us define the shape functions for the linear axisymmetric triangular elements with coordinates (R i , Z i ), (R j , Z j ), (R k , Z k ) and the area |A e |. So the linear shape functions are where Using (21)- (23), the first local stiffness matrix is and A e is the area of the triangular element. The second stiffness component comes from the first boundary integral in (20), so we have where L ij denotes the length of the side of the triangular element connecting vertices i and j. For evaluation of the last component of the stiffness matrix we apply the Gauss-Legendre quadrature and use seven Gaussian points (2n − 1 = 12 and n = 6.5), so we have where the subscript p represents the evaluation of the corresponding function at the pth Gaussian point, w p are Gaussian weights, R p = (L 1 ) p R i + (L 2 ) p R j + (1 − L 1 − L 2 ) p R k , and |J| is the Jacobian defined by Now we evaluate two load vectors which are the second integral in (20) and the integral on the RHS of (19). So, we have The other one is also evaluated by the Gauss-Legendre approach, hence The global matrix Equation using (25)-(30) is K{T} = {F}.

Description of the Parametric Study
As an example, temperature results for the parallel-plate using the finite difference model were obtained for a fixed distance from the center to the wall y of 1.9 mm and the width z of 200 cm. Figure 3a presents the temperature distribution along (z) and the half of the distance (y) between parallel-plates for a case with the wall temperature (T wall ) of 433.15 K and the inlet temperature (T inlet ) of 403.15 K. As can be noticed, the temperature of the flow gets more evenly distributed along the y axis as z increases (See Figure 3b).  Also the analyses were conducted with the finite element model considering the 12,482 triangular elements to determine the temperature distribution along the length of the tube. As an example, Figure 4a shows the temperature distribution along the length (Z) and the radius (R) of the tube for a case with the wall temperature (T wall ) of 433.15 K and inlet temperature (T inlet ) of 403.15 K. The locations along the R and Z were normalized by the radius (2.2 mm) and the length (800 mm) of the tube respectively. As seen in Figure 4b, the temperature of the polypropylene melt flow gets more uniform along the direction of the R as the Z increases, and then finally reaches a steady state distribution.
The required length of parallel-plate or circular tube is defined as the length of the channel in the direction of the flow from the entrance such that the temperature of the flow reaches the steady state level. To quantify the required length of the tube (L req ) to reach the steady state temperature distribution, the following criteria were used: where i = 0, 1, . . . is associated with the location along R and l = 0, 1, . . . is associated with the location along Z in finite element solution. When the criterion in (31) is satisfied, the steady state location along Z will be l + 1 and the required length is identified. For the example shown in Figure 4, the l is identified as 25 and the corresponding L req is determined as 260 mm.

Results and Discussion
Following the same procedure, a parametric study was conducted for wall temperatures, inlet temperatures, radius of the tube, and pressure drop to determine the L req . Figures 5 and 6 show the parametric results for the flow temperature in the parallel-plate and circular tube channels. In general, the L req increases as the radius of the tube increases since the larger radius can accommodate more fluid volume. As seen in Figures 5a and 6a, when the inlet temperature increases and approaches the wall temperature, the L req becomes smaller. The explanation of this is that for the higher inlet temperature, it requires less distance for the flow temperature to reach the steady state. On the other hand (see Figures 5b and 6b), as the wall temperature increases, the L req becomes larger since it needs more distance for the flow temperature to reach the uniform distribution. Also note, the L req increases when the pressure drop increases since larger pressure means higher flow speed (see Figures 5c and 6c). The radius of the tube has the highest effect on the required length of the tube while other parameters have minor effects, especially for the pressure drop and the wall temperature.
It is important to note the manifolds with decreasing depth or also constant depth, which is the subject of our study. For instance, Karkri and Jarny [33] analyzed the temperature profile in the parallel-plate with constant depth using the conjugate gradient method that shows the results similar to our solutions. Another paper by Andreozzi et al. presented the temperature profile results of the fluid flows in parallel-plate systems and the authors confirmed them by the Ansys-Fluent output [34]. Their results on temperature profiles of the fluid showed similar behavior as in our results ( Figure 3). As it was mentioned before, the finite element model was suitable to solve the temperature profile in a circular tube. In comparison to the solutions obtained by Wei and Luo [18], our temperature results were consistent (Figure 4).
The parametric analyses would be helpful in the designing of the extrusion die channels as these parametric study results may contribute in proper selection (for instance, radius of the tube) or controlling the parameters (i.e., wall or inlet temperature) and suggest the optimal length of the channels. Eventually, the optimal length derived from the study would assure an efficient production in which correspondingly the flow temperature is evenly distributed.

Conclusions
Casson rheology models were investigated for the velocity and temperature profiles, particularly, of polypropylene melt flows in extrusion die channels. Finite difference and finite element methods were used to solve for temperature distributions of the flow in the channels. The numerical model provides estimates of locations in the channels where the steady state temperature distribution is reached in Casson fluid melts. The results in the parametric study indicate that the required parallel-plate and tube lengths are highly affected by the width of the channels while other parameters (wall temperature, inlet temperature and pressure drop) have minor effects. These results can be used as a base for constructing the design of the die and increase the production efficiency of PP sheets and PP films. In addition, the numerical approach can be applied for other materials other than polypropylene, for instance, the materials that follow the Casson rheology model. Further research can be conducted by applying other rheology models such as Cross or Carreau-Yasuda for better understanding the nature of the PP flows.
T wall wall temperature u velocity of the fluid flow y axis coordinate perpendicular to z in parallel-plate y c ,y crit critical length from the center to the critical point along y in parallel-plate y max distance from the center to the wall in parallel-plate z, Z dimensional, and non-dimensional axial coordinates in the direction of the flow F (e) ,{F} local and global load vectors K (e) ,K local and global stiffness matrices W vector of three linear shape functions Superscripts and subscripts (e), e evaluation of the entity referring to the particular element p evaluation of the corresponding function at the particular point T matrix transpose J Jacobian matrix