Modeling melt conveying and power consumption of conveying elements in co-rotating twin-screw extruders

1 Introduction

The co-rotating fully intermeshing twin-screw extruder is one of the most commonly used extruder types, as it provides excellent mixing in combination with sufficient pressure build-up. The conveying elements are responsible for most of the pressure generation and forward conveying, and are therefore an integral part of the screw configuration of this extruder type. The conveying behavior of this screw element directly influences the back-pressure length (fully filled region), which is needed to generate the pressure at the screw tip caused by the die. The power-consumption behavior, in contrast, is directly linked to the melt temperature at the screw tip. Conveying elements are commonly modeled by a flat-plate model as described by White and Kim (2010), who detailed the development of twin-screw modeling over the years. Note that the flat-plate model is only an approximation and does not capture the entire complex 3D geometry with all clearances and the intersection area. The intersection area, for example, has been previously represented by a shift of the screw channel (Szydlowski and White 1987) and by an additional stationary section with adapted geometry (Potente et al. 1999). Szydlowski and White (1987) illustrated the effect of including the intersection area in their approximation by directly comparing their results to those of Denson and Hwang (1980), who ignored the intersection area. In our previous work (Stritzinger et al. 2023a), we additionally investigated in detail and without any simplifications the influences of the clearance sizes between barrel and screw and between the two screws. For this purpose, we performed CFD simulations without geometry simplifications and varied the clearance sizes over an extensive range of pitches and diameter ratios. Conveying elements are based on the self-wiping Erdmenger profile, which was first discussed by Erdmenger (1964) and later described mathematically by Booy (1978). To minimize the number of influencing parameters, we transferred the complex 3D geometry of conveying elements into dimensionless space and showed that this geometry is fully described by four dimensionless geometry parameters (Stritzinger et al. 2023a): the dimensionless diameter ratio Π _D, the dimensionless pitch Π _T, the dimensionless screw clearance Π _δ, and the dimensionless nip clearance Π _s, which are the input parameters for our models. For the transition of the geometry and process parameters the Buckingham Π-theorem according to Durst (2008) was used with barrel diameter D, screw speed N, and shear viscosity η as basic quantities. To this end, the screw pitch T, the screw-screw clearance s, and the screw-barrel clearance δ are normalized with the barrel diameter D. Moreover, only for the outer screw diameter D _a and the core screw diameter D _k the well know diameter ratio Π _D is used. Furthermore, the definition of the dimensionless flow volume Π _V was adapted to mirror the volume of a cylinder and not the volume of a sphere, thus the length of the conveying element L was included. The flow volume along V is the free volume in the extruder between the two conveying elements and the barrel and therefore is a function of the other geometry parameters. According to these guidelines, the following definitions of the dimensionless geometry parameters were derived:

To describe the flow in conveying elements, also dimensionless process parameters are needed. Consequently, the volume flow rate V ˙ , the axial pressure gradient ∂p/∂z, the mechanical drive power P, and the viscous dissipation Q ˙ Diss are transformed into the dimensionless space with the same basic quantities. This gives the following definitions for the dimensionless process parameters:

Due to the linear relationship of both the dimensionless pressure gradient Π p ′ and the dimensionless mechanical power Π _Pow to the dimensionless volume flow-rate Π V ˙ for Newtonian fluids, the target parameters of our models are profile parameters that, according to the theory of similarity, describe these linear functions, as reported by Kohlgrüber et al. (2020a).

To achieve an increased compatibility with models of commonly used kneading blocks (Stritzinger et al. 2023b), the following adapted profile parameters were chosen as target parameters: dimensionless drag-flow capacity A ₁, dimensionless element conductance A ₃, dimensionless turning point B ₂, and dimensionless turbine parameter B ₃. With these dimensionless profile parameters, the conveying and power-consumption behaviors can be calculated, respectively, by

The dimensionless element conductance A ₃ is defined as the ratio between the dimensionless drag-flow capacity A ₁ and the dimensionless dam-up pressure A ₂. Consequently, it can be expressed as the reciprocal slope of the linear function between the dimensionless pressure gradient Π p ′ and the dimensionless volume flow-rate Π V ˙ . The dimensionless turbine parameter B ₃ quantifies the ratio between the dimensionless turning point B ₂ and the dimensionless turbine point B ₁. It is the slope of the linear function that relates the dimensionless mechanical power Π _Pow to the dimensionless volume flow-rate Π V ˙ . All of the mentioned profile parameters (A ₁. A ₂, A ₃, B ₁, B ₂, and B ₃) are only depending on the geometry of the conveying element.

2 Dataset

As a basis for symbolic regression, a vast number of CFD simulations were carried out with Ansys Polyflow 2019 R2 (ANSYS Polyflow 2019) to determine the conveying and power-consumption parameters for various dimensionless geometry parameters. For the simulations, the complex 3D geometry of the conveying elements, as illustrated in Figure 1, was considered, and we used the same solver settings and simulation set-up as in our previous work (Stritzinger et al. 2023a). For each geometry, two simulations were carried out: one to determine the drag-flow capacity directly, with two moving screws with a screw speed of N = 100  rpm and zero pressure gradient between the periodic in- and outlet; and one to directly determine the element conductance, with stationary screws and a defined volume flow rate of V ˙ = 10 − 4   m 3   s − 1 through the periodic in- and outlet. The dimensionless power parameters were calculated by linear superposition, which we described in great detail in our previous work on kneading blocks (Stritzinger et al. 2022).

Figure 1:

Fluid domain used for CFD simulations including all clearances and the intersection region. The colors indicate the boundary conditions: Green – periodic inlet (Δp = 0 / m ˙  = fixed), red – barrel (stationary wall), yellow – screw 1 (rotational velocity = 100  rpm / rotational velocity = 0), and blue – screw 2 (rotational velocity = 100  rpm / rotational velocity = 0).

To generate a training dataset, we carried out a comprehensive subfactorial design study in the parameter space given in Table 1. By using wide value ranges for each dimensionless geometry parameter, we sought to include the dimensions of all commonly used conveying elements.

We selected 772 design points (∼50 %) from the entire full factorial set for our parameter study to limit the computational time required. To ensure that the entire parameter space was sampled and no knowledge was lost, we included the center points of the parameter space, the edges and surfaces and the corner points (as illustrated in Figure 2) and chose the remaining design points randomly in the hypercube to avoid bias.

Figure 2:

Schematic drawing of the parameter space, including the edge center points (yellow), the surface center points (orange), and the corner points (red).

Further, a validation dataset was generated using 100 randomly chosen design points within the parameter space that were not already included in the training dataset. This dataset was used to test our models’ interpolation capability (i.e., its ability to predict previously unseen conveying elements within the parameter space).

The influence of the geometry parameters on the profile parameters according to the simulation dataset was already discussed in our previous work (Stritzinger et al. 2023a) and can be summed up with Table 2. The clearances play only a minor role when compared to the dimensionless pitch and diameter ratio. To highlight the relationship between the influencing and target parameters and to sum up our previously published work, the Pearson correlation coefficient R and the coefficient of determination R² are evaluated and depicted in Figure 3.

Figure 3:

Strength of linear correlation between geometry and profile parameters.

Accordingly, to R and R² the linear relationship between the influencing and target parameters is evaluated and ranked. The coefficient of determination is only above 0.3 for five combinations of influencing and target parameters and thus symbolic regression is a good option to encompass all of the information contained in the dataset in easy-to-use algebraic functions which outperform linear regression. Each of the target parameters can be roughly estimated by at least one influencing parameter with the following linear functions:

3 Symbolic regression

Using the training dataset, symbolic regression models were trained that generated easy-to-use analytical formulas from the knowledge gained from the CFD simulations. For this purpose, the software Heuristic Lab (Wagner et al. 2014) was used to simultaneously optimize model accuracy and complexity.

3.1 Algorithm settings

The multi-objective NSGA-II algorithm was applied and model complexity evaluated in terms of tree size (see Figure 4) and in terms of accuracy by Pearson’s R². The maximum tree size was set to 100, and a population size of 500, a maximum generation number of 500 and a mutation rate of 25 % were defined for the algorithm. The model grammar included addition, subtraction, multiplication, division, and root, quadratic and cubic functions for all models. The default tree grammar was expanded with exponential and logarithmic functions if the target accuracies of a mean relative error of less than 1 % for the conveying parameters and less than 3 % for the power parameters were not met. For each target parameter at least 10 symbolic regression models were trained, then the best mean-relative-error model was chosen.

Figure 4:

Schematic of a symbolic regression model structure and example for two different complexities.

4 Comparison to state-of-the-art approaches

To further validate the accuracy of our models, we compared these to measurements (Düphans et al. 2024), CFD simulations (Düphans et al. 2024; Kohlgrüber et al. 2020b) and models of the conveying behavior published by other research groups (Potente et al. 1990, 1994). For this comparison with state-of-the-art analytical, numerical and experimental approaches, we used the three geometries given by Düphans et al. (2024) and listed in Table 4, and generated the characteristic screw lines for each model and simulation result. As an analytical method, the flat-plate model by Potente et al. (1990, 1994) defined in Equations (20) and (21) with the geometry relationships provided in the Appendix was used:

The two numerical approaches we chose for comparison were the CFD simulations carried out by Düphans et al. (2024) using Ansys Fluent 2020 (ANSYS Inc., Canonsburg, PA, USA) and the simulation results presented by Kohlgrüber et al. (2020b) for typical two-lobe extruder elements, where no clear information was given on the clearance sizes or simulation method or program used for the graphs presented. Düphans et al. (2024) used the moving-wall boundary condition for the screw surfaces and included a die at the end of the conveying elements with three different diameters. With this simulation set-up, they avoided transient simulations with dynamic meshing and determined the conveying and power parameters for the three conveying element geometries. They additionally designed a test rig and used silicone oil to experimentally determine the conveying and power parameters. However, they did not measure the radial temperature in their experiments, which we would have needed for proper evaluation of the power parameters, as these are closely linked to dissipation. We thus compared only their conveying parameters to those obtained from our models. The values from the literature and from our hybrid models are plotted in Figure 5.

Figure 5:

Comparison of literature values (from measurements, simulations, and analytical modeling) to the results of our prediction models for three conveying elements with different screw pitches.

Especially, the critical region near the dam-up pressure was predicted most accurately by our regression models. Conveying elements are commonly deployed within an operating window close to dam-up pressure, which makes this region the most significant in the graphs. Due to the slightly more pronounced influence of channel curvature, the flat-plate model provided poorer predictions for lower pitches. The simulation results of the research group in Dortmund (Düphans et al. 2024) deviated markedly from those presented by Kohlgrüber et al. (2020b). Presumably, this is due to the moving-wall condition in Fluent, which models tangential rotational motion accurately only if the geometry is a surface of revolution (e.g., a cylinder) (ANSYS Fluent 2020). The difference between our results and the simulation results by Kohlgrüber et al. (2020b) are probably due to the higher clearance sizes used by Kohlgrüber et al. (2020b), as pointed out in Chapter 4 of his book, where he presents simulation results of conveying elements with screw-barrel clearances of 0.5 mm and 1 mm. According to the results of our parameter study, greater screw-barrel clearance leads to a decreasing drag-flow capacity, as already discussed in our previous work (Stritzinger et al. 2023a).

5 Application of models for non-isothermal shear thinning material flow

To highlight the application of our models in real life polymer processing, the example of the calculation of a back-pressure length in a real-life compounding application is presented in this section. A simple set-up for the production of a polymer blend is used for this purpose and presented in Figure 6.

Figure 6:

Schematic depiction of maximum back-pressure length l _max to avoid flooding of the degassing opening of a co-rotating twin screw extruder.

Our models can be used to estimate the difference between the maximum back-pressure length l _max and the actual back-pressure length l _full and therefore can help avoid flooding of the degassing opening. The profile parameters A ₁, A ₃, B ₂, and B ₃ can be calculated using the geometry of the conveying elements between the end of the degassing zone to the screw tip and combined with process conditions like throughput m ˙ , screw speed N, measured melt temperature T _m, and pressure p _die at the screw tip. To calculate the length of the fully-filled region in front of the die l _full the throughput m ˙ must first be transformed into a volume flow-rate V ˙ and then converted into the dimensionless volume flow-rate Π V ˙ according to Equation (6). For this transformation the melt density of the polymer needs to be calculated according to the Tait equation (Osswald and Hernández-Ortiz 2006):

The dimensionless pressure gradient can be calculated using the dimensionless drag-flow capacity A ₁ and the dimensionless element conductance A ₃, as shown in Equation (10). In order to transform the pressure gradient back into dimensional space and calculate the back-pressure length, it is also necessary to calculate the shear viscosity of the polymers for the melt temperature and a representative shear rate. For the representative shear rate our previously presented thermal representative shear rate γ ˙ rep , therm (Stritzinger et al. 2023b) can be used in combination with the approximated Arrhenius model (Rauwendaal 2014) and the Carreau model (Osswald and Hernández-Ortiz 2006) for the viscosity:

Applying Equation (7) to our use-case the following formular for the actual back-pressure length is reached:

This is just one simple example of how the models can be used for real-life compounding applications. An additional very important use case is the modeling of backward-conveying elements at the end of melting sections, which guarantee fully-filled kneading blocks in the melting zones and are therefore crucial for the performance of the melting zone and the sealing of the degassing zone from the hopper.

6 Conclusions

We generated symbolic regression models that describe the conveying and power-consumption behaviors of fully intermeshing co-rotating twin-screw conveying elements. Our models exploit the accuracy of the 1,544 CFD simulations on which they are based, and they can easily be applied to all commonly used conveying elements. A direct comparison to state-of-the-art models and experiments proves the power of our approach and underlines the improvement in twin-screw extrusion modeling. Our prediction models could now be used in a variety of scenarios, and deploying their accurate inclusion of all clearances, curvatures, and intersection areas will stimulate novel applications.