Numerical Simulation of Polymer Solutions in a Single-Screw Extruder

Abstract

Single-screw extruders are the most common equipment used for polymer extrusion. The study of the hydrodynamics of a polymer melts flow in the extruder channel is the basis for modeling and understanding the extrusion process. In general form, the extruder includes a straight section with a screw installed in it. In this study, the three-dimensional mathematical modeling of the polymer solutions flow in the metering zone of a single-screw extruder is performed. The influences of the screw geometry (L/D₂ = 1…3) on the flow structure and the pressure drop are analyzed under a speed rotation up to 60 rpm. Aqueous solutions of 0.5% polyacrylamide (0.5% PAA) and 1.5% sodium salt of carboxymethyl cellulose (1.5% CMC) are considered as the working fluid.

1. Introduction

The single-screw extruders are the classic equipment in the processing of polymers [1]. Films, sheets, pipes, and various profiles are produced by the extrusion method. Simulation of the extrusion process allows us to study the influence of many parameters (such as geometrical, rheological) on the structure of the polymer melt and predict the final product [2]. It is very difficult to observe the processes occurring in the extruder channel, therefore one of the possible ways to control and determine the rational parameters of the entire extrusion process is mathematical modeling.

In [3,4], it is noted that currently, it is difficult to obtain an exact analytical solution to even a simplified mathematical model of the extrusion process (one-dimensional non-Newtonian flow). Analytical solutions of the problem of the Newtonian fluid isothermal flow in the channel of a single-screw extruder are known, which differ in the way they solve the system of differential equations [5,6]. Several works are devoted to semi-analytical solutions of a similar problem, where the rheological behavior of the polymer is described by a power law [7,8,9].

Experimental studies of polymers (polymer melts and solutions) revealed the nonlinear viscoelastic behavior [10,11,12,13,14]. Simulation of polymer flow (polymer solutions and melts) including viscoelastic effects is a difficult task. In this regard, in most published works, elastic effects are not taken into account and polymeric media are considered as shear-thinning fluids. The power law model is the most frequently used model for the rheological description of polymer solutions and melts, for example, [2,15,16].

Despite the complexity of the problem being solved, the application of numerical methods is an urgent task. Numerical methods make it possible to predict the structure of the polymer melt flow, including the presence of normal stresses, adjust the equipment operating parameters, and reduce the cost of experimental studies.

When developing mathematical models, the most commonly used assumptions are those that reduce the original problem to a two-dimensional formulation. For example, the assumption that the helical screw channel is unwound and considered as a flat rectangular channel [4,17,18,19,20,21,22]. Thus, plane-parallel fluid flow in a channel of infinite width is considered instead of direct 3D modeling. Also, the principle of reversed motion is often used (tube rotation is considered at the fixed screw) [4,17,18,19,20,21,22]. In all works, polymers are considered as shear-thinning fluid [4,17,18,19,20,21,22].

The main functional zones of the extruder are feeding (solids conveying), melting, and metering. The metering zone is the simplest for modeling in a single-screw extruder, since the flow laws of viscous fluids can be applied to the polymer melt flow [23]. In this zone, the polymer is completely molten [24]. Besides, it is the metering zone that determines the extruder performance [23].

Work [4] describes a simple numerical method for solving the problem of the polymer melt flow in the metering zone of the single-screw extruder. The 2D mathematical model takes into account the presence of circulating motion, non-Newtonian behavior of the processed material, and heating of the polymer due to viscous dissipation.

In [3,18,19,21,22,25], numerical solutions obtained in a 2D and 3D formulation are approximated by symbolic regression models. The heuristic approach for predicting the main characteristics of the polymer melts flow in the extruder channel is based on genetic programming. Speaking about the evolutionary computation in general, it should be noted that they, like any method using an element of randomness, do not guarantee the detection of a global extremum of the objective function (or optimal solution) for a certain time. Their main advantage is that they allow finding a “better” solution to a difficult task in less time than other methods.

In contrast with the standard flat-plate model of the unwound screw channel (2D), this hydrodynamic problem is solved in a three-dimensional formulation in [17,26]. Authors of the papers consider mathematical modeling of the flow in the metering zone of a plasticizing extruder. The objective of the study was the melting of the polymer composition based on polyethylene.

Work [27] presents the results of three-dimensional numerical simulation of the polymer melt flow in a channel with a screw (extruder). In this case, the polymer melt was considered as the Newtonian fluid, which does not correspond to the real rheological behavior.

In the model presented in [28], the continuous penalty finite element scheme and the power-law model were used to solve the basic equations of continuity and momentum in a three-dimensional Cartesian coordinate system.

Work [2] presents the three-dimensional modeling of the wood-plastic composites melt flow in the metering section of the single-screw extruder. The paper notes that the extrusion of these polymers is significantly different from the extrusion of plastics. Therefore, it is important to consider that the rheological behavior of different polymers can vary significantly.

Some results of three-dimensional modeling of heat and mass transfer under laminar flow of the low-density polyethylene in a single-screw extruder is presented [29]. The Carreau-Yasuda model, which includes five parameters, was used to describe the rheology behavior of the polymer melt. However, the drawback of this work is the lack of analysis of mesh convergence; numerical results were obtained using a coarse mesh.

Work [24] discusses the study of secondary flows (Moffatt eddies) arising in the metering zone of the extruder. Open-source computational fluid dynamics package OpenFoam was used for numerical studies. The rheological properties of shear-thinning fluid were described with the Carreau model.

The article [16] presents a three-dimensional model of a single-screw extruder. The partial rotational fraction has been implemented with the finite element. However, as noted by the authors, the extension of the developed method in viscoelastic flow simulations is not possible.

The aim of this work is to study the structure of the polymer solutions flow in a single screw extruder using 3D modeling. The 0.5% polyacrylamide (0.5% PAA) and 1.5% sodium salt of carboxymethyl cellulose (1.5% CMC) are used as working fluids. This study is part of an ongoing study whose goal is to develop a mathematical model to describe the flow of viscoelastic fluids in a single-screw extruder.

2. Simulation Methodology

2.1. Geometric Representation, Governing Equations, and Boundary Conditions

This work considers a stationary laminar three-dimensional flow of shear-thinning fluid in the metering zone of a single-screw extruder. The part of the extruder under investigation is presented in the form of a channel with a screw twisted counterclockwise (Figure 1). Mathematical modeling was performed under the assumption that gravity force is negligible and the fluid is incompressible [21]. In this regard, the distributions of hydrodynamic characteristics are identical in areas #1 and #2. Therefore, further, area #1 is considered as the area under investigation (Figure 2). In this paper, the screw rotation is simulated by the rotation of the outer pipe, this method is widely used for similar problems in the two-dimensional setting [19].

The computational area includes a straight section of length $L_0$ and a channel section with a screw of length $L_1$. It should be noted that the outer wall rotation is organized in section $L_1$. Here L is the pitch of the screw channel, i.e., the length of the channel part that corresponds to the rotation of the screw flights on 360 degrees, h is the flight thickness, $r_1$, $r_2$ are inner and outer radii of the channels.

The governing equations in tensor form are as follows:

(1)$\rho \left(\nabla v\cdot v\right)=\nabla \cdot \left(-\rho I+\eta \left(\nabla v+{\left(\nabla v\right)}^T\right)\right)$

(2)$\nabla \cdot v=0$

where p—pressure, Pa; I—unit tensor; η—effective viscosity, Pa·s; $\rho$—fluid density, kg/m³; $v=\left\{v_r,v_{\phi },v_z\right\}$—velocity vector.

Boundary conditions:

(1) At the inlet, a developed velocity profile in the coaxial channel is used [30]

(3)$v_z=2v_a\frac{\left(1-{\left(r/r_2\right)}^2\right)\log \left(r_1/r_2\right)-\left(1-{\left(r_1/r_2\right)}^2\right)\cdot \log \left(r/r_2\right)}{\left(1+{\left(r_1/r_2\right)}^2\log \left(r_1/r_2\right)+\left(1-{\left(r_1/r_2\right)}^2\right)\right)},v_r=0,v_{\phi }=0$

where $v_a$—average velocity, defined to be the volumetric flow rate divided by the cross-sectional area of flow;

(2) The no-slip conditions of the fluid are assumed on the channel walls of straight section (L₀) and on the surfaces of the screw flight and screw shaft ($r=r_1$) [16,19,21,29,31];

(3) On the tube wall ($r=r_2$):

(4)$v_{\phi }=\left(-1\right)\cdot r_2\cdot 2\cdot \pi \cdot N$

where $N$—a number of full rotations per sec, therefore $N^{*}=N\cdot 60$—a number of full rotations per min. Let us note that because of the adopted orientation of the coordinate system (Figure 2), the $v_{\phi }$ component is taken with a negative sign to implement the correct direction of rotation of the outer pipe.

The numerical solution of the problem is implemented in the commercial package COMSOL Multiphysics^®. It is known that, in this package, the basic Equations (1) and (2) include only the pressure gradient while the pressure function itself is absent. Thus, the pressure is determined accurately up to a constant value. Based on the above, the condition that pressure is equal to zero is specified as the boundary condition at the channel outlet [22].

2.2. Rheological Properties of Polymer Solutions

This study examined aqueous solutions of polyacrylamide and sodium carboxymethyl cellulose. PAA is used very widely because of its gelling property. For example, it is used in the process of wastewater treatment, in paper production, extraction of mineral resources, oil production, agriculture, medicine, construction, and other industries. Polyacrylamide is one of the best polymers used as a drag reduction agent and vortex inhibitor because of its flexibility and long chain [32].

CMC is currently the most widely used biopolymer in the world made of cellulose. Owing to its availability, it is also in great demand in various industries, for example, in oil and gas production, mining-and-chemical industry, production of detergents, as a stabilizer, thickener, filler in the food industry, pharmaceutical industry, etc. Besides, it is used as a binder in the production of coal dust briquettes, pencils, and chalks, as well as in granulation of seeds [33].

Polymer solutions have a nonlinear viscoelastic flow pattern. The paper does not consider viscoelastic effects because of the complexity of the problem being solved. This assumption is quite common [34].

Cross model was used to describe the rheological behavior of these fluids [35,36]:

(5)$\eta ={\eta }_{\infty }+\frac{{\eta }_{{}_0}-{\eta }_{\infty }}{1+{\left(\lambda \dot{\gamma }\right)}^n}$

where ${\eta }_0$, ${\eta }_{\infty }$ Newtonian plateaus (asymptotic values of dynamic viscosity at $\dot{\gamma }\to 0$ and $\dot{\gamma }\to \infty$ respectively); ${\lambda }_c$—specific time; n—coefficient determining the slope of the viscosity curve; $\dot{\gamma }$—shear rate, 1/s. It should be mentioned that the Cross model was also used in the following studies [17,32,37,38].

The parameters of the Cross model for 1.5% CMC and 0.5% PAA are tabulated in Table 1. The experimental data obtained from [39] served as the basis for obtaining the parameters of the Cross model for 0.5% PAA.

Figure 3 shows the curves of viscosity versus shear rate, where the experimental data were taken from [39].

2.3. Mesh Parameters

For numerical simulation, it is important to determine the mesh convergence of the obtained results by studying the convergence according to the computational mesh. Hence, we performed additional numerical studies and determined the optimal mesh sizes. The studies were performed by the example of an isothermal flow of 0.5% PAA in the channel with a screw. It was previously noted that we are considering one of the areas bounded by a screw (Figure 1b and Figure 2). Figure 4 presents the computational area. The geometric parameters of the computational area: $L/D_2$ = 2, $D_2$ = 0.039 m, $D_1=0.026$ m, $L_0=1.5\cdot D_2$, $L_1=5.5\cdot D_2$, h = 0.004 m.

It is known that the axial component of the velocity vector is less sensitive to the mesh parameters than the radial or tangential one. Therefore, the influence of the mesh parameters on the numerical results was estimated by comparing the profiles of the radial component of the velocity vector around the circumference in cross-sections $A_{45}$ and $A_{360}$ (dashed line in Figure 4b). Here $A_i$ (i = 45…360) is the cross-section spaced apart from $A_0$ at a distance corresponding to the rotation of the screw by i degrees, $r_3=\left(r_1+r_2\right)/2$, β is the current angle value, α is the sector size. The problem was numerically solved by using a server with two Intel Xeon Gold 6126 processors (24 cores in total, 512 Gb RAM).

Figure 5 presents the profiles of the dimensionless radial component of the velocity vector for the three partitions of the computational area. The partition of the computational area into 132,413 (mesh #1), 543,657 (mesh #2), 1,051,349 (mesh #3) elements was considered. It can be seen that the profiles of the radial component of the velocity vector coincide when using mesh #2 and mesh #3. Even though the profile $v_r/v_a$ for mesh #1 practically coincides with the results obtained on a finer mesh, there is a small region in the vicinity of the right fin in which $v_r/v_a$ is nonzero, which contradicts the no-slip boundary conditions of the fluid.

The performed analysis of mesh convergence showed that mesh #2 parameters provide sufficient numerical results that are independent of further partitioning of the computational area. Table 2 shows the parameters of the used mesh. Since we considered laminar flows of a viscous fluid at low Reynolds numbers (Re < 11) then the channel length $L_1=5.5\cdot D_2$ with $L/D_2$ = 2 pitch of twist of screw is sufficient to exclude output effects from the channel.

2.4. Validation

The problem of a 0.2% CMC laminar flow in an annular channel was considered for the validation of the numerical results and the solution method implemented using the commercial package COMSOL Multiphysics^®. The geometric parameters of the computational area correspond to $D_2$ = 0.067 m, $D_1$ = 0.032 m, $v_a$ = 0.69 m/s [40]. Figure 6 presents the comparison results of the profile of the dimensionless axial velocity. In contrast to the results of [40], the obtained velocity profile is consistent with the boundary conditions on the channel walls, namely, the no-slip conditions of fluid, i.e., $v_z$ = 0.

In addition, the comparison between numerical results (profile of axial velocity) and distribution (3) is performed for water flow in a coaxial channel with following parameters (Figure 7): 8D—the length of coaxial channel; $D_2$ = 0.039 m, $D_1$ = 0.026 m, $V_a$ = 0.05 m/s. In addition, the profiles of the axial velocity for 0.5% PAA and 1.5% CMC are shown in the same figure. For 0.5% PAA, the velocity profile is more flattened than the parabolic profile related to water. This is explained by the fact that in the near-wall region, the viscosity of 0.5% PAA is much lower than in the central part of the channel.

3. Results and Discussion

This section presents the numerical studies results of the 0.5% PAA and 1.5% CMC flows in a channel with a screw. The geometric characteristics of the channel (Figure 1) are as follows: the outer diameter of the pipe $D_2$ = 0.039 m, the inner diameter $D_1$ = 0.026 m, the length of the straight section $L_0=1.5\cdot D_2$. The following values of the pitch of twist of the screw fins $L/D_2$ = 1, 2, and 3 were considered. The screw is twisted counterclockwise, and the length of the channel with the screw depends on the pitch of twist and varies from 3$D_2$ to 7.5$D_2$ depending on the pitch of twist ($L/D_2$). The main calculations were performed for an average velocity of the fluid of 0.05 m/s.

Figure 8 shows the lines of equal values of the dimensionless axial (Figure 8a) and radial (Figure 8b) velocity, as well as viscosity (Figure 8c) for 0.5% PAA flow in the channel with a screw. Designations $A_{45}$, $A_{90}$, $A_{360}$ are cross-sections of the channel according to Figure 4. Figure 8c shows the dimensionless viscosity, which is the ratio of the actual viscosity to the average viscosity ${\eta }_{inlet}$ calculated in the inlet cross-section of the channel. For 0.5% PAA ${\eta }_{inlet}$ = 0.0509 Pa·s. As can be seen from the figure, when polymer solution flows into the curved section, provided by the presence of a screw, the velocity and viscosity profiles are rearranged. Areas with the maximum axial velocity appear at the inlet to the screw, then as the fluid flows along the channel, the velocity profile becomes more flatted. At the average velocity of 0.05 m/s, it was revealed that a steady-state flow regime is achieved in the section up to cross-section $A_{360}$. Besides, the steady-state flow regime of 0.5% PAA is characterized by the presence of a radial velocity (up to 11%) comparable to the average velocity. This indicates the presence of radial fluid flows. At a steady-state flow, the flatted velocity profile is accompanied by the presence of an area with very high viscosity (Figure 8c) (a double exceedance of viscosity with respect to the average viscosity at the channel inlet).

Profiles of the dimensionless axial velocity in various cross-sections of the channel (from A₀ to A₄₅₀) along a circle with a radius $r=r_3$ (see Figure 4) were plotted to determine the length of the initial hydrodynamic section in the channel with a screw (Figure 9). Results are presented for the case when outer pipe (barrel) is rotated with speed N* = 0, 30, 60 rpm. It was found that for 1.5% CMC, the hydrodynamic entrance length is bounded by the cross-section A₁₈₀, while for 0.5% PAA, it is bounded by the cross-section $A_{360}$. This is explained by the fact that at the same average velocity of two fluid, the corresponding Reynolds numbers differ by two orders of magnitude, namely Re (0.5% PAA) = 10.64, Re (1.5% CMC) = 0.87. It should be noted that the Reynolds number was calculated using the average viscosity (${\eta }_{inlet}$) in the inlet section of the channel and the hydraulic diameter $D_h=4S/p$, where $S$ is the cross-sectional area of the channel, $p$ is the wetted perimeter.

Figure 10 shows the axial velocity profiles for the steady-state flow of the studied fluids for various rotation velocities of the outer pipe ($r=r_3$, see Figure 4). For both fluids, an increase in the outer tube rotation velocity leads to a slight deformation of the velocity profile, which is consistent with the results of the work [16]. It can be seen that the greatest deformation is observed for 0.5% PAA since its apparent viscosity is less than that of 1.5% CMC. It should be noted that for 0.5% PAA the axial velocity profile has the inflection point near the left screw flight caused by the rotation velocity of the outer pipe. Thus, the obtained results indicate that the outer pipe rotation makes the velocity profiles more flatted.

Figure 11 shows the axial velocity profiles for fully developed flow in the channels with the different pitch of twist of the screw ($L/D_2$). Results are presented for equal average velocity (Figure 11a) and equal Reynolds number (Figure 11b). The corresponding profiles were plotted along a circle with a radius $r=r_3$. It can be seen (Figure 11a) that with an increase of the twist intensity (decrease of the $L/D_2$ parameter), the velocity profile for 1.5% CMC becomes more elongated. Another behavior is observed for the flow of 0.5% PAA. With a decrease of the $L/D_2$ parameter, the velocity profile takes on a more flatted form. So in the near-wall region the velocity of 0.5% PAA increases compared with 1.5% CMC. The observed increase of velocity in the near-wall region is explained by the presence of higher velocity gradients, which contribute to a decrease in viscosity. A similar result was obtained for the same Reynolds numbers (Figure 11b), in this case, the displacement of the axial velocity profile to the left wall is pronounced for 0.5% PAA, where the direction of displacement coincides with the direction of rotation of the outer pipe.

Figure 12 shows the dependence of the pressure drop on the number of revolutions of the outer pipe and average velocity at a fully developed flow of 0.5% PAA and 1.5% CMC in the metering zone of a single-screw extruder (distance between A₃₆₀ and A₄₅₀ cross-sections). It was found that with an increase in the $L/D_2$ parameter (with an increase of the twist intensity), a decrease in the pressure drop due to an increase in the number of revolutions is observed for both fluids (Figure 12a). An increase in the average velocity leads to an increase in the pressure drop (Figure 12b). Processing the numerical results shown in Figure 12 revealed that the dependences of the pressure drop on the number of revolutions and average velocity could be approximated by a second-order polynomial with accuracy of approximation $R^2$ = 0.999. For example, for 0.5% PAA, the following relationships were obtained:

$\begin{array}{l}L/D_2=1:dP/dz=-5787.6\cdot {\left(N*\right)}^2+19708\cdot N*+436.06(R^2=0.999)\\ L/D_2=2:dP/dz=-75389\cdot {\left(N*\right)}^2+15599\cdot N*-74.412(R^2=0.999)\\ L/D_2=3:dP/dz=-54978\cdot {\left(N*\right)}^2+11439\cdot N*-88.396(R^2=0.999)\end{array}$

It can be seen from the figure (Figure 12a,b) that in both cases the highest values of the pressure drop correspond to 1.5% CMC (blue lines), which is due to its higher viscosity compared to that of 0.5% PAA (see Figure 3).

4. Conclusions

Mathematical modeling of the steady-state flow of shear-thinning fluid in the metering zone of a single-screw extruder was performed. The different geometries of screw were chosen to investigate the flow structure of polymer solutions in extruder with rotating outer tube. The flows of aqueous polymer solutions of 0.5% PAA and 1.5% CMC were considered. Numerical results showed that for 1.5% CMC, the hydrodynamic entrance length in the extruder (speed rotating up to 60 cycle/min) does not exceed the distance limited by a 180-degree rotation of the screw, and for 0.5% PAA, it does not exceed the distance limited by a 360-degree rotation of the screw. The obtained numerical results showed that, for the aqueous polymer solutions flow in the metering zone of the extruder, the profile of the axial velocity is more flatted for 0.5% PAA than for 1.5% CMC because of its rheological behavior and the presence of rotational motion in the extruder. The velocity profile becomes more flatted as the rotational speed becomes greater. It is observed that an increase in the rotation velocity leads to a decrease in the pressure drop, which facilitates the fluid flow along the extruder length. It was found that the dependence of the pressure drop on the number of revolutions as well as on the average velocity can be fitted by a second-order polynomial function with accuracy $R^2$ = 0.999.

Author Contributions

A.K., R.G., and J.K. conceived and developed the paper; R.G. performed numerical simulations and contributed to manuscript preparation; J.K. wrote the review and contributed to manuscript preparation; A.K. provided guidance throughout the development of this paper.

Funding

The research was conducted with financial support by RSF, project number 19-11-00220.

Conflicts of Interest

The authors declare no conflicts of interest.

Figures and Tables

Figure 1. Transversal (a) and longitudinal (b) channel section.

Figure 2. Computational area.

Figure 3. The dependence of viscosity on shear rate.

Figure 4. Computational area (a) and area of the working fluid flow in cross-section A0 (b).

Figure 5. (a) Radial velocity distribution (va = 0.05 m/s) and (b) fragment of the partition of the computational area into triangular elements (mesh #2).

Figure 6. Profile of the dimensionless axial velocity in the channel cross-section: 1—data [40]; 2—COMSOL Multiphysics®.

Figure 7. Profile of the dimensionless axial velocity in the coaxial channel.

Figure 8. Lines of equal values of dimensionless axial (a) and radial (b) velocity, (c) viscosity (L/D2 = 2, N* = 0).

Figure 9. Profiles of the dimensionless axial velocity in various cross-sections of the channel (Va = 0.05 m/s): 1—A0, 2—A45, 3—A90, 4—A135, 5—A180, 6—A360, and 7—A450 are dashed line in black.

Figure 10. Profiles of the dimensionless axial velocity for different the rotational speed N*, (L/D2 = 2, va = 0.05 m/s).

Figure 11. Profiles of the dimensionless axial velocity in the cross-section A360 for 0.5% polyacrylamide (PAA) (red) and 1.5% carboxymethyl cellulose (CMC) (blue): 1 and 3—L/D2 = 1 (dashed line); 2 and 4—L/D2 = 3 (solid line).

Figure 12. Pressure drop for fully developed flow in extruder for 0.5% PAA (red) and 1.5% CMC (blue).

Table 1

Cross model—characteristic material parameters at 294 K.

Name	Concentration	η0 (Pa·s)	η∞ (Pa·s)	λc (s)	n (-)	ρ (kg/m3)
1.5% CMC [35]	1500 ppm	1.97	0.02	0.172	0.64	1025
0.5% PAA	500 ppm	1.08	0.0023	5	0.8	999

Table 2

Parameters of mesh #2.

Parameters	Subdomain	Boundary
Maximum element size	0.0013	0.0011
Type of elements	Triangle (regular)	Triangle (advancing front)
Minimum element quality	0.3469	0.3469