*1.1. Analysis of Flow in Metering Channels*

Early theoretical analyses that investigated melt conveying in single-screw extruders dealt with one- and two-dimensional flows of temperature-independent Newtonian

**Citation:** Marschik, C.; Roland, W.; Dörner, M.; Steinbichler, G.; Schöppner, V. Leakage-Flow Models for Screw Extruders. *Polymers* **2021**, *13*, 1919. https://doi.org/10.3390/ polym13121919

Academic Editors: Salah Aldin Faroughi, Luís L. Ferrás, Alexandre M. Afonso and Célio Bruno Pinto Fernandes

Received: 4 May 2021 Accepted: 7 June 2021 Published: 9 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

fluids. These problems have exact analytical solutions for the drag and pressure flows in the cross- and down-channel directions, respectively. The first model of screw viscosity pumps was published anonymously [8] and later extended by Rowell and Finlayson [9]. Similarly, Carley et al. [10] proposed a simplified flow theory for screw extruders. Mohr et al. [11,12] investigated the characteristics of the transverse flow in the screw channel. Both the complexity and the accuracy of the analysis increase when the non-Newtonian flow behavior of polymer melts is included. It is well known that polymer melts are shear-thinning fluids whose viscosity decreases with increasing shear rate. Pseudo-plastic behavior complicates the mathematical model such that the governing flow equations must be solved numerically, and exact closed-form analytical solutions are no longer possible. The viscosity being dependent on the shear rate, the drag and pressure flows are coupled. For multidimensional flows, complexity is increased further by the combined effect of shear in the down- and cross-channel directions. To gain insights into how a material's shear-thinning nature affects melt conveying in single-screw extruders, several authors presented numerical solutions for one- and two-dimensional flows of power-law fluids in metering channels. Rotem and Shinnar [13] obtained numerical results for a one-dimensional flow under isothermal conditions. Including the effect of transverse flow, Griffith [14], Zamodits and Pearson [15], and Karwe and Jaluria [16] presented numerical solutions for shear-thinning flows in infinitely-wide screw channels.

With the development of more advanced computers, a trend emerged towards a numerical analysis of three-dimensional flows in metering channels. Including non-Newtonian and non-isothermal effects, a number of studies provided detailed insights into the recirculating extruder channel flow. Spalding et al. [17] applied the finite-element method to examine the flow in a helical metering section. Ghoreishy et al. [18] presented a non-isothermal analysis of elastomer melt flow in unwound screw channels. Polychronopoulos and Vlachopoulos [19] analyzed the development of Moffat Eddies in the flight root corners of the screw channel. Vachagina et al. [20] investigated the non-isothermal conveying behavior of a Giesekus fluid in a helical screw section.

To avoid numerical procedures, various approximate solutions for shear-thinning flows were developed. Booy [21] employed the classical Newtonian pumping model to derive effective viscosity values and approximate shear-thinning flow behavior. Rauwendaal [22] proposed correction factors to the drag- and pressure flows in the tra-ditional model to include the shear-rate-dependent flow behavior of polymer melts. Similarly, Spalding et al. [23] introduced a correction factor for the drag flow.

Rather than presenting corrections to the traditional Newtonian theory, we pro-posed an alternative approach for power-law fluids. Using a hybrid modeling approach that combines analytical, numerical, and data-based modeling techniques, we approximated a large number of numerical solutions of scaled flow equations obtained for two- and three-dimensional formulations. Assuming two-dimensional flows of power-law fluids in infinitely-wide screw channels, Pachner et al. [24] and Roland et al. [25,26] proposed regression models for the prediction of throughput and viscous dissipation. To include the influence of the flight flanks, we [27–30] extended the theories to three-dimensional flows. These widely applicable models increase prediction accuracy by including the combined effects of shear-thinning flow behavior, transverse flow, and screw flights.

#### *1.2. Analysis of Leakage Flow*

The influence of leakage flow on the pumping capability of extruder screws has been investigated since the earliest theories for the metering zone. Rowell and Finlayson [9] approximated the annular clearance between flight land and barrel surface by two parallel plates and assumed the leakage flow to be a pressure flow through an infinitely wide slit. Gore and McKelvey [31] extended the analysis by adding the effect of flight clearance on the drag flow. Mohr and Mallouk [12] presented a more sophisticated theory by considering the pressure flow caused by the cross-channel pressure gradient in the screw channel and the effect of flight clearance on the drag flow. A detailed review of the early Newtonian theories

was provided by Tadmor and Klein [1]. For extruder screws with "normal" clearances, the velocity profile in the flight clearance was shown to be determined mainly by the drag-flow component and influenced only marginally by the transverse pressure gradient.

Taking the shear-thinning flow behavior of polymer melts into account, Rauwendaal and Housz [32] used a numerical approach to readdress the flow characteristics in the flight clearance. For large flight clearances and small power-law indices, a significant influence of the transverse pressure flow was identified. Further, the total power consumption of the screw was shown to be affected considerably by the power consumption in the clearance. To account for the effect of leakage on throughput, corrections were made to (i) the transverse flow in the screw channel and (ii) the net flow rate. Focusing on pressuregenerating metering zones, their analysis omitted overridden functional zones that are subject to a negative axial pressure gradient.

Including non-isothermal effects, Meyer et al. [33] numerically evaluated the temperature development in the flight clearance. Taking a pure drag flow into account, the results indicated that the thermal development length is generally smaller than the available gap length. When predicting the velocity and temperature distribution at the exit of the flight clearance, the flow can therefore be considered fully developed. To investigate the effect of leakage flow on the temperature distribution in the screw channel, Pittman and Rashid [34] and Rauwendaal [35] carried out numerical analyses of two- and three-dimensional power-law-model-based flows.

#### *1.3. Research Approach*

This work presents a new modeling approach to including leakage effects in the analysis of melt-conveying zones. Rather than directly correcting the drag and pressure flows, we developed new analytical regressions to predict the flow rate and dissipation locally in the leakage gap. In the construction of the leakage-flow models, particular attention was paid to their usability in our screw-simulation routine for predicting the conveying behavior of melt-conveying zones. For detailed information, the reader is referred to the following articles and PhD-theses: Aigner [36], Marschik et al. [37–39], Roland et al. [40], and Luger [41].

The main idea of the simulation approach is based on network theory: To reduce the complexity of the flow situation, the screw channel is subdivided into very small interconnected segments of constant cross-section and physical parameters. These sections are represented by network elements that consist of a source and a resistance connected in parallel. Using melt-conveying models from the literature, the network elements describe the flow locally in the screw channel. Analogously to electrical circuits, the network elements are connected via nodal points to form an equivalent flow network, whose flow rates and pressure differences are evaluated by means of nodal analysis. In contrast to timeconsuming numerical simulations based on computational fluid dynamics, the method iteratively solves a linearized set of network equations and is therefore significantly faster.

To illustrate the objective of this research, Figure 1 shows a flow network of a conventional metering zone, which consists of down- and cross-channel elements. Due to the changes in channel depth in the transverse direction, the latter is initialized with three network elements in series: (a) one from the channel center to the pushing flight, (b) one over the flight clearance, and (c) one from the trailing flight to the channel center one turn behind. The extrusion literature provides numerous theories for predicting the characteristics of the down-channel flow. Accurate models for the analysis of leakage effects, in contrast, remain elusive.

A flow chart of our work is given in Figure 2. First, we described mathematically the flow under consideration, converted the model into a dimensionless form, and carried out a detailed analysis of the governing equations. The dimensionless model was then solved numerically to evaluate flow rate and dissipation for a large number of physically independent setups. To remove the need for numerical methods, we derived approximate equations for flow rate and dissipation, using symbolic regression based on genetic programming. The structure of the regression models and their parameters were optimized iteratively until we obtained satisfying prediction accuracies. The new regressions were designed to be able to describe the properties of the network elements positioned in the flight gap.

**Figure 1.** The development of leakage-flow models to predict flow rate and dissipation in the leakage gap of a discretized screw channel.

**Figure 2.** Schematic work-flow chart of the hybrid modeling approach, including analytical, numerical, and data-based modeling.

#### **2. Model Development**

We used a hybrid modeling approach to derive approximate equations for the flow rate and the dissipation in the flight clearance. The method incorporates (i) analytical, (ii) numerical, and (iii) data-based modeling techniques. A detailed description of the main characteristics of our hybrid modeling approach was provided in [42].

#### *2.1. Analytical Modeling*

#### 2.1.1. Problem Definition

We applied the flat-channel approximation to simplify the helical reference frame of the screw. The screw channel is unwound and considered as a rectangular flow channel covered by an infinite flat plate. The top view of the unwound screw channel is shown in Figure 3. The simplified channel representation is based on a Cartesian coordinate system: *x* is the cross-channel direction and *z* the down-channel direction. Avoiding cylindrical coordinates, the effect of channel curvature is ignored. The validity of the flat-plate model was discussed critically by Sun and Rauwendaal [43] and Roland et al. [44] for Newtonian and non-Newtonian polymer melts, respectively. Both studies confirmed the correctness of the approximation for small channel-depth-to-diameter ratios. This geometrical precondition is particularly fulfilled in the flight clearance. For conventional extruder screws with standard clearances, the ratio is typically in the range of *δ*/*D<sup>b</sup>* = 0.001. As a result, the error introduced by unwinding the screw channel at the radial position of the barrel is small. The influence of channel curvature can hence be ignored. The following parameters were used to describe the geometry of the screw: *D<sup>b</sup>* is the screw diameter, *t* the screw pitch, *e* the flight width, *w* the channel width, and *ϕ<sup>b</sup>* the pitch angle.

**Figure 3.** The top view of the unwound screw channel.

Drawing on the traditional pumping model, we reversed the kinematics of the process; that is, the screw is stationary, and the barrel moves at circumferential speed v*<sup>b</sup>* :

$$\mathbf{v}\_b = D\_b \,\,\pi \,\,\mathbf{N}\_\prime \tag{1}$$

where *N* is the screw speed. The barrel velocity can be decomposed into a cross- and a down-channel component, v*b*,*<sup>x</sup>* and v*b*,*<sup>z</sup>* :

$$\mathbf{v}\_{b,\mathbf{x}} = \mathbf{v}\_b \sin(\varphi\_b); \qquad \qquad \mathbf{v}\_{b,\mathbf{z}} = \mathbf{v}\_b \cos(\varphi\_b); \tag{2}$$

$$
\tan(\varphi\_b) = \frac{t}{D\_b \ \pi}.\tag{3}
$$

Unlike previous studies that investigated the effect of the flight clearance [9,12,31,32], we locally analyzed the flow between flight land and barrel surface. Figure 4 illustrates the flat-plate model of the screw channel and the flight clearance.

**Figure 4.** Flat-plate approximation of the unwound screw channel (**a**) and the flight clearance (**b**). In the representation of the flight clearance, the flow channel was rotated clockwise by 90◦ : *x* is the direction across and *z* the direction along the clearance.

The space through which the leakage flow occurs can be considered infinitely wide, resulting in a two-dimensional flow. Note that the coordinate system is the same in both cases. For convenience, however, the flow channel of the flight clearance was rotated clockwise by 90◦ : (i) *x* is the direction across the clearance, and (ii) *z* is the direction along the clearance. The channel gap between barrel surface and flight land is given by the clearance *δ*.

#### 2.1.2. Governing Equations

In the first step, we derived the governing equations to mathematically describe the flow through the flight clearance. To this end, we made the following assumptions: (i) the flow is independent of time, fully developed, and isothermal, (ii) the fluid is incompressible and wall-adhering, and (iii) gravitational forces are ignored. A detailed discussion on the validity of these commonly-applied assumptions was provided by Tadmor and Klein [1]. The velocity field was reduced to a two-dimensional flow:

$$\mathbf{v} = \begin{pmatrix} \mathbf{v}\_x(y) \\ \mathbf{0} \\ \mathbf{v}\_z(y) \end{pmatrix}. \tag{4}$$

Assuming a fully developed flow, the velocity components are functions of the channelheight coordinate *y* only; that is, the two-dimensional formulation is not capable of predicting any flow rate variations in the *x*- and *z*-directions. Meyer et al. [33] confirmed the validity of the assumption for both Newtonian and power-law fluids. Due to the low velocities and high viscosities of extruder flows, the Reynolds number is usually small (Re1), and the flow can be considered laminar [28,30]. Viscous forces being dominant over inertial forces, we reduced the problem to Stokes flow. By omitting time derivates, we ended up with the following simplified momentum equations:

$$\frac{\partial p}{\partial x} = \frac{\partial \tau\_{yx}}{\partial y}. \tag{5}$$

$$\frac{\partial p}{\partial z} = \frac{\partial \pi\_{yz}}{\partial y}. \tag{6}$$

In our isothermal theory, the momentum equations are uncoupled from the energy equation, which allowed us to evaluate viscous dissipation from the velocity field subject to the following boundary conditions:

$$\mathbf{v}\_x(y=0) = \mathbf{0}.\tag{7}$$

$$\mathbf{v}\_x(y=\delta) = \mathbf{v}\_{b,\mathbf{x}}.\tag{7}$$

$$\mathbf{v}\_z(\mathcal{y}=\mathbf{0})=\mathbf{0}.\tag{8}$$

$$\mathbf{v}\_z(\mathcal{y}=\delta)=\mathbf{v}\_{b,z}.\tag{8}$$

The rheological behavior of the polymer melt was expressed by a nonlinear constitutive equation that relates the stress tensor *τ* to the rate-of-deformation tensor *D*, which is given by the symmetric part of the velocity-gradient tensor *L*:

$$
\boldsymbol{\pi} = \mathbf{2} \,\eta \,\mathbf{D}.\tag{9}
$$

$$\mathbf{D} = \frac{1}{2} \left( \mathbf{L} + \mathbf{L}^{\mathsf{T}} \right). \qquad \qquad \qquad \mathbf{L} = \nabla \mathbf{v}. \tag{10}$$

For most polymer melts, the viscosity is strongly dependent on the velocity gradients in the flow field, rendering them shear thinning. To describe the shear-rate dependency of the viscosity, we applied a power-law model, where *K* is the consistency and *n* the power-law index. The magnitude of shear rate . *γ* is related to the second invariant of the rate-of-deformation tensor:

$$
\eta = K|\dot{\gamma}|^{n-1}.\tag{11}
$$

$$\left|\dot{\gamma}\right| = \sqrt{2(\mathbf{D}:\mathbf{D})} = \left[\left(\frac{\partial \mathbf{v}\_x}{\partial y}\right)^2 + \left(\frac{\partial \mathbf{v}\_z}{\partial y}\right)^2\right]^{\frac{1}{2}}.\tag{12}$$

With these definitions, the shear stresses were obtained from:

$$\tau\_{y\chi} = \eta \frac{\partial \mathbf{v}\_{\chi}}{\partial y} = K \left[ \left( \frac{\partial \mathbf{v}\_{\chi}}{\partial y} \right)^{2} + \left( \frac{\partial \mathbf{v}\_{z}}{\partial y} \right)^{2} \right]^{\frac{n-1}{2}} \frac{\partial \mathbf{v}\_{\chi}}{\partial y}. \tag{13}$$

$$\pi\_{yz} = \eta \frac{\partial \mathbf{v}\_z}{\partial y} = K \left[ \left( \frac{\partial \mathbf{v}\_x}{\partial y} \right)^2 + \left( \frac{\partial \mathbf{v}\_z}{\partial y} \right)^2 \right]^{\frac{n-1}{2}} \frac{\partial \mathbf{v}\_z}{\partial y}. \tag{14}$$

Finally, the momentum equations were rewritten as:

$$\frac{\partial p}{\partial x} = \frac{\partial}{\partial y} \left\{ K \left[ \left( \frac{\partial \mathbf{v}\_x}{\partial y} \right)^2 + \left( \frac{\partial \mathbf{v}\_z}{\partial y} \right)^2 \right]^{\frac{n-1}{2}} \frac{\partial \mathbf{v}\_x}{\partial y} \right\}. \tag{15}$$

$$\frac{\partial p}{\partial z} = \frac{\partial}{\partial y} \left\{ K \left[ \left( \frac{\partial \mathbf{v}\_x}{\partial y} \right)^2 + \left( \frac{\partial \mathbf{v}\_z}{\partial y} \right)^2 \right]^{\frac{n-1}{2}} \frac{\partial \mathbf{v}\_z}{\partial y} \right\}. \tag{16}$$

These coupled nonlinear partial differential equations are the governing equations of our model. In combination with the boundary conditions in Equations (7) and (8), the equation system describes the velocity field of the two-dimensional flow. Since, generally, valid closed-form analytical solutions are unknown, the shooting method was used to calculate the velocity field and the volume flow rate . *V* per unit width and total viscous dissipation per unit area (. *qdiss* = *τ* : *L*):

$$
\dot{V} = \int\_0^\delta \mathbf{v}\_x(y) dy. \tag{17}
$$

$$\dot{q}\_{\rm diss} = \eta \left[ \left( \frac{\partial \mathbf{v}\_x}{\partial y} \right)^2 + \left( \frac{\partial \mathbf{v}\_z}{\partial y} \right)^2 \right]. \tag{18}$$

$$
\dot{Q}\_{\rm diss} = \int\_0^\delta \dot{q}\_{\rm diss}(y) dy. \tag{19}
$$

#### 2.1.3. Theory of Similarity

Before the flow equations were solved numerically, the model was converted into a dimensionless form, using the theory of similarity, for the following reasons [26–30]: (i) Two systems described by the same dimensionless quantities are similar; (ii) varying any one of the independent input parameters changes the physical conditions of the flow; (iii) solutions obtained for a specific set of dimensionless input parameters apply to all dimensional variations that yield the same dimensionless input parameters.

We introduced the following dimensionless spatial coordinate and fluid velocities:

$$\zeta = \frac{y}{\delta}. \qquad \qquad \qquad v\_{\chi} = \frac{\mathbf{v}\_{\chi}}{\mathbf{v}\_{b,\chi}} \qquad \qquad \qquad v\_{z} = \frac{\mathbf{v}\_{z}}{\mathbf{v}\_{b,\chi}}. \tag{20}$$

These definitions were used to derive dimensionless parameters for shear rate and viscosity in the flight clearance:

$$\eta^\*\_{\delta} = \frac{\eta \cdot \delta^{n-1}}{K \operatorname{v}^{n-1}\_{b,\mathbf{x}}} = \left| \dot{\gamma}^\* \right|\_{\delta}^{n-1} = \left[ \left( \frac{\partial v\_{\mathbf{x}}}{\partial \xi} \right)^2 + \left( \frac{\partial v\_{\mathbf{z}}}{\partial \xi} \right)^2 \right]^{\frac{n-1}{2}}.\tag{21}$$

Similarly, the momentum equations were transformed into dimensionless form:

$$6\operatorname{II}\_{p,\mathbf{x}}^{\delta} = \frac{\partial}{\partial\xi} \left\{ \left[ \left( \frac{\partial v\_{\mathbf{x}}}{\partial\xi} \right)^{2} + \left( \frac{\partial v\_{z}}{\partial\xi} \right)^{2} \right]^{\frac{n-1}{2}} \frac{\partial v\_{\mathbf{x}}}{\partial\xi} \right\};\tag{22}$$

$$6\operatorname{II}\_{p,z}^{\delta} = \frac{\partial}{\partial \overline{\xi}} \left\{ \left[ \left( \frac{\partial v\_x}{\partial \overline{\xi}} \right)^2 + \left( \frac{\partial v\_z}{\partial \overline{\xi}} \right)^2 \right]^{\frac{n-1}{2}} \frac{\partial v\_z}{\partial \overline{\xi}} \right\} \tag{23}$$

where the parameters Π*<sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* and Π*<sup>δ</sup> <sup>p</sup>*,*<sup>z</sup>* are dimensionless pressure gradients:

$$\Pi\_{p,i}^{\delta} = \frac{\frac{\partial p}{\partial \mathbf{i}} \, \delta^{n+1}}{6 \, \mathbf{K} \, \mathbf{v}\_{b,\mathbf{x}}^{n}} \, \mathbf{y} \qquad \mathbf{i} = \mathbf{x} \, \mathbf{z}. \tag{24}$$

In addition, the boundary conditions were rewritten to:

$$v\_{\mathfrak{x}}(\mathfrak{x}=\mathbf{0})=\mathbf{0}.\qquad\qquad\qquad v\_{\mathfrak{x}}(\mathfrak{x}=\mathbf{1})=\mathbf{1}.\tag{25}$$

$$v\_z(\xi=0) = 0. \qquad \qquad v\_z(\xi=1) = \tan(90^\circ - \varphi\_b). \tag{26}$$

Solving the dimensionless momentum Equations (22) and (23) in combination with the boundary conditions (25) and (26) required an additional mathematical constraint to be defined. Previous theories that modeled two-dimensional flows of polymer melts in infinitely wide screw channels considered the cross-channel net flow to be zero [22,26]. Physically, this means that leakage flow across the screw flights was ignored. Since this assumption is not reasonable in the analysis presented here, we instead assumed the pressure gradient in the *z*-direction to be zero. In our local formulation, this implies that the corresponding flow component is governed exclusively by the rotation of the screw:

$$
\frac{\partial p}{\partial z} = 0 \quad \rightarrow \quad \Pi^\delta\_{p,z} = 0. \tag{27}
$$

Our assumption is based on the results of a similar study [26], in which the authors examined the flow of polymer melts in infinitely wide screw channels. For a variety of screw designs and processing conditions, the dimensionless down-channel pressure gradient along the screw channel *Πp*,*<sup>z</sup>* was shown to be within the range of (−1.0; 1.0). Rather than using the channel depth *h*, our theory requires the flight clearance *δ* to be predefined for the calculation of the dimensionless pressure gradient in the *z*-direction. Since the latter is significantly smaller, the dimensionless down-channel pressure gradient along the flight clearance *Π<sup>δ</sup> p*,*z* converges to zero.

A comparison of typical parameter values of *Πp*,*<sup>z</sup>* and *Π<sup>δ</sup> p*,*z* is given in Table 1. The results indicate the different orders of magnitude of the parameters.


**Table 1.** Dimensionless pressure gradients in the screw channel *Πp*,*<sup>z</sup>* and in the flight clearance *Π<sup>δ</sup> p*,*z*

Finally, we derived dimensionless parameters for the volume flow rate *Π<sup>δ</sup> V* and the total viscous dissipation *Π<sup>δ</sup> Q* :

$$\Pi\_V^\pounds = \frac{2\not\dot{V}}{\delta\operatorname{v}\_{b,\text{x}}} = 2\int\_0^1 v\_{\text{x}}(\xi) \, d\xi. \tag{28}$$

.

$$II\_Q^\delta = \frac{\dot{Q}\_{\rm diss} \,\delta^n}{\text{K }\mathbf{v}\_{b,\mathbf{x}}^{n+1}} = \int\_0^1 \left[ \left(\frac{\partial v\_\mathbf{x}}{\partial \boldsymbol{\xi}}\right)^2 + \left(\frac{\partial v\_z}{\partial \boldsymbol{\xi}}\right)^2 \right]^\frac{n+1}{2} d\boldsymbol{\xi}.\tag{29}$$

Our dimensionless model has three independent input parameters that completely describe the physics of the flow: (i) the screw-pitch ratio *t*/*D<sup>b</sup>* , indicated by *tan*(*ϕb*), (ii) the power-law index *n*, and (iii) a dimensionless pressure gradient across the flight clearance *Π<sup>δ</sup> p*,*x* . The first is part of the boundary conditions (25), while the second and third are included in the momentum Equations (22) and (23).

2.1.4. Set-Up of Parametric Study

We created four sets of physically independent design points by varying the input parameters *t*/*D<sup>b</sup>* , *n*, and *Π<sup>δ</sup> p*,*x* (Tables 2–4). For all sets, the screw-pitch ratio and the powerlaw index were varied within the ranges (0.5–2.47) and (0.2–1.0), respectively. These ranges include almost all extruder screws and polymer melts in industrial use. The variation of the dimensionless pressure gradient across the flight clearance *Π<sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* was adjusted case by case. To determine the industrially relevant parameter ranges, we carried out screw calculations for various types of extruder screws, using our network-based routine [1–40].

**Table 2.** The ranges of variation for *t*/*D<sup>b</sup>* , *n*, and *Π<sup>δ</sup> p*,*x* for Data Set 1.



**Table 3.** The ranges of variation for *t*/*D<sup>b</sup>* and *n* for Data Sets 2 to 4.

**Table 4.** The ranges of variation for *Π<sup>δ</sup> p*,*x* for Data Sets 2 to 4.


For conventional extruder screws with normal clearances (*δ* ≈ 0.001 *D<sup>b</sup>* ), the dimensionless pressure gradient is relatively small. Its magnitude, however, increases if the extruder screw is constructed with an undercut flight. This situation can be found, for example, in barrier screws. High dimensionless pressure gradients arise in the context of high-performance screws, such as energy-transfer screws, in which screw flights are, by design, undercut to promote cross-channel mixing. In addition, the direction of leakage flow depends on the conveying characteristics of the screw. For pressure-generating functional zones, leakage flow typically reduces the throughput. For strongly overridden functional zones, in contrast, polymer melt passing through the clearance is forced towards the screw tip to increase the net flow rate of the processing machine. This effect can be increased if the flight clearance is undercut.

For Data Set 1, the dimensionless pressure gradient was varied within the range (−1.0; 1.0). Since pressure development also depends on the shear-thinning nature of the polymer melt, the minima, maxima, and increments for Data Sets 2–4 were based on the power-law index. In total, we created 11,781 design points. Finally, Data Sets 1–4 were merged, and design points with multiple occurrences were deleted, which yielded 9231 physically independent setups.

#### *2.2. Numerical Modeling*

#### 2.2.1. Numerical Solution Procedure

In the next step, we numerically evaluated the target variables (dimensionless volume flow rate *Π<sup>δ</sup> V* and dimensionless dissipation *Π<sup>δ</sup> Q* ) of our model for all 9,231 physically independent setups, using the shooting method. To this end, we calculated the velocity field in the flight clearance by solving the governing equations of our dimensionless model. A detailed description of the numerical solution procedure was given in [26]. Transforming the boundary value into an initial-value problem, we derived explicit forms of the dimensionless momentum equations:

$$\frac{\partial v\_{\chi}}{\partial \tilde{\xi}} = \frac{1}{\eta\_{\delta}^{\*}} \Big( \delta \, \Pi\_{p,\mathbf{x}}^{\delta} \tilde{\xi} + \mathbb{C}\_{1} \Big);\tag{30}$$

$$\frac{\partial v\_z}{\partial \xi} = \frac{\mathbf{C}\_2}{\eta\_\delta^\*},\tag{31}$$

where the dimensionless viscosity was rewritten as:

$$\eta\_{\delta}^{\*} = \left[ \left( \theta \, \Pi\_{p,x}^{\delta} \xi + \mathcal{C}\_1 \right)^2 + \mathcal{C}\_2^2 \right]^{\frac{n-1}{2n}}.\tag{32}$$

The initial estimates of the integration constants *C*<sup>1</sup> and *C*<sup>2</sup> were taken from the Newtonian solution. Applying the Simpsons rule yielded the following equations for the velocity profiles:

$$v\_{\mathfrak{x}}(\mathfrak{x}) = \underbrace{v\_{\mathfrak{x}}(\mathfrak{x} = \mathbf{0})}\_{0} + \int\_{0}^{1} \frac{\partial v\_{\mathfrak{x}}}{\partial \mathfrak{x}} d\mathfrak{x}.\tag{33}$$

$$v\_z(\mathfrak{f}) = \underbrace{v\_z(\mathfrak{f}=\mathbf{0})}\_{0} + \int\_0^1 \frac{\partial v\_z}{\partial \mathfrak{f}} d\mathfrak{f}.\tag{34}$$

To iteratively solve the unknowns, we used a Newton-Raphson scheme:

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - \mathbf{J}(\mathbf{x}\_n)^{-1} [\mathbf{f}(\mathbf{x}\_n) - \mathbf{f}(\mathbf{x})],\tag{35}$$

where *x* is the vector of unknowns, *J* the Jacobian matrix, and *f* is the vector of boundary conditions:

$$
\begin{pmatrix} \mathbb{C}\_{1,n+1} \\ \mathbb{C}\_{2,n+1} \end{pmatrix} = \begin{pmatrix} \mathbb{C}\_{1,n} \\ \mathbb{C}\_{2,n} \end{pmatrix} - \begin{pmatrix} \frac{\partial v\_{\ge 1}}{\partial \mathbb{C}\_1} & \frac{\partial v\_{\ge 1}}{\partial \mathbb{C}\_2} \\ \frac{\partial v\_{\ge 1}}{\partial \mathbb{C}\_1} & \frac{\partial v\_{\ge 1}}{\partial \mathbb{C}\_2} \end{pmatrix}^{-1} \begin{bmatrix} \begin{pmatrix} v\_{\ge 1,n} \\ v\_{\ge 1,n} \end{pmatrix} - \begin{pmatrix} \tan(\rho\_b) \\ 1 \end{pmatrix} \end{pmatrix}. \tag{36}
$$

The velocity boundary conditions at the barrel surface will not be met unless the initial values are perfect. The converged solutions for the velocity profiles were then used to determine the dimensionless target variables *Π<sup>δ</sup> V* and *Π<sup>δ</sup> Q* .

For all calculations, the dimensionless channel height was divided into 1000 equidistant segments. A solution was considered converged if the difference in dimensionless volume flow rate *Π<sup>δ</sup> V* between two iterations was smaller than 10-8. Previous analyses have shown that these settings are sufficient to obtain mesh-independent results for our target variables [26].

#### 2.2.2. Numerical Results

Our parametric design study encompassing 9,231 independent setups provided numerical solutions for the dimensionless volume flow rate *Π<sup>δ</sup> V* and the dissipation *Π<sup>δ</sup> Q* in the flight clearance as functions of the dimensionless input parameters *t*/*D<sup>b</sup>* , *n*, and *Π<sup>δ</sup> p*,*x* .

Figure 5 shows the dimensionless volume flow rate as a function of the dimensionless pressure gradient across the flight clearance for various power-law indices. For all setups, the curves are symmetrical about the point of pure drag flow (*Π<sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* = 0); that is, an equidistant increase in the pressure gradient (positive or negative) affects the magnitude of the dimensionless volume flow rate equally.

The power-law index is a measure of the shear-thinning behavior of the polymer melt: the lower the power-law index, the more shear-thinning is the fluid. Assuming a Newtonian fluid, the widely-known linear behavior is evident: For *Π<sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* = 0 (pure drag flow), the curve satisfies *Π<sup>δ</sup> <sup>V</sup>* = 1, while for *<sup>Π</sup><sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* = 1, the zero-throughput condition is fulfilled *Π<sup>δ</sup> <sup>V</sup>* = 0.

**Figure 5.** The volume flow rate *Π<sup>δ</sup> V* as a function of pressure gradient *Π<sup>δ</sup> p*,*x* : Influence of power-law index for *t*/*D<sup>b</sup>* = 0.5 (**a**) and *t*/*D<sup>b</sup>* = 1.0 (**b**).

Generally, the volume flow rates become negative if a critical pressure gradient is reached. From a mathematical viewpoint, this means that the direction of flow changes to the negative *x*-direction. Physically, it implies that the pressure flow caused by the pressure build-up across the clearance exceeds the drag flow, which yields a negative net throughput. According to our model definition, this behavior is subject to strongly overridden melt-conveying zones. Positive flow rates, in contrast, are found in pressuregenerating metering zones, in which the leakage flow reduces the net throughput. With decreasing power-law index, the critical pressure gradient shifts to lower values.

The curves become increasingly nonlinear and pressure-sensitive with decreasing power-law index. Two effects are evident: (i) For a given dimensionless volume flow rate, the more the dimensionless pressure gradient (positive and negative) increases, the less shear-thinning the fluid. (ii) For highly shear-thinning polymer melts, small variations in the pressure gradient can lead to pronounced variations in the volume flow rate.

The influence of the screw-pitch ratio on the dimensionless volume flow rate is less pronounced (Figure 6). For positive pressure gradients, the target variable increases with increasing screw-pitch ratio, while the opposite behavior is evident for negative pressure gradients. This effect is caused by the influence of the flow along the flight clearance (in the *z*-direction) on the deformation rates, which becomes more pronounced the lower the screw-pitch ratio.

For markedly positive or negative dimensionless pressure gradients (Figure 6b), the effect of the screw-pitch ratio decreases significantly since the flow is governed mainly by the pressure gradient across the flight clearance.

Figure 7 illustrates the influence of the power-law index on the dimensionless dissipation for a square-pitched screw with *t*/*D<sup>b</sup>* = 1. Viscous dissipation is mainly responsible for the temperature development in the channel. Due to inner friction, mechanical energy is converted into heat, causing a rise in melt temperature.

**Figure 6.** The volume flow rate *Π<sup>δ</sup> V* as a function of pressure gradient *Π<sup>δ</sup> p*,*x* : Influence of screw-pitch ratio for *n* = 0.2. The scaling of the diagrams was adjusted to better visualize the influence for smaller pressure gradients (**a**) and larger pressure gradients (**b**).

**Figure 7.** Dimensionless dissipation *Π<sup>δ</sup> Q* as a function of the dimensionless pressure gradient *Π<sup>δ</sup> p*,*x* (**a**) and as a function of the dimensionless volume flow rate *Π<sup>δ</sup> V* . (**b**) The influence of the power-law index for a square-pitched screw with *t*/*D<sup>b</sup>* = **1.0**.

The dimensionless dissipation can be plotted as a function of either the dimensionless pressure gradient (Figure 7a) or the dimensionless volume flow rate (Figure 7b). For both representations, the curves are again symmetrical about the point of pure drag flow (*Π<sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* = 0 or *Π<sup>δ</sup> <sup>V</sup>* = 1), where the target variable reaches a minimum. In general, the dimensionless dissipation increases if the pressure flow contributes to the flow characteristics; that is, the higher the dimensionless pressure gradient, the more pronounced is the frictional heat generation. Similarly, dissipation becomes highly dependent on the power-law index for strongly pressure-generating or pressure-consuming flows.

Three effects are observed: (i) For moderate pressure gradients, the more dimensionless dissipation decreases, the more shear-thinning the polymer melt. (ii) For higher magnitudes, in contrast, the more frictional heat generation decreases, the more Newtonian the fluid. (iii) For constant dimensionless volume flow rates, viscous heating increases with the increasing power-law index.

Figure 8 shows the influence of the screw-pitch ratio on the dimensionless dissipation for a polymer melt with power-law index n = 0.2. For both constant dimensionless pressure gradients (Figure 8a) and constant dimensionless flow rates (Figure 8b), viscous

dissipation increases with decreasing screw-pitch ratio. This result is again caused by the effect of transverse flow in the leakage gap. The influence of the screw-pitch ratio on dimensionless dissipation vanishes almost completely in strongly pressure-generating and pressure-consuming flows.

**Figure 8.** Dimensionless dissipation *Π<sup>δ</sup> Q* as a function of the dimensionless pressure gradient *Π<sup>δ</sup> p*,*x* (**a**) and as a function of the dimensionless volume flow rate *Π<sup>δ</sup> V* . (**b**) The influence of the screw-pitch ratio for a polymer melt with power-law index *n* = **0.2**.

To represent the diverse characteristics of the leakage flow, we considered a wide range of dimensionless pressure gradients. Especially for highly shear-thinning polymer melts with low power-law indices, our extended dataset caused a significant nonlinear increase in the target variables, yielding values higher than *<sup>Π</sup><sup>δ</sup> V* <sup>&</sup>gt; <sup>10</sup><sup>5</sup> and *Π<sup>δ</sup> Q* <sup>&</sup>gt; <sup>10</sup><sup>6</sup> , as illustrated in Figure 9. Recently, we have shown that the parameters are limited to *Π<sup>V</sup>*  < 40 and *Π<sup>Q</sup>* <sup>&</sup>lt; <sup>140</sup> when analyzing the flow in metering channels [26]. This comparison illustrates the increased complexity of the following symbolic regression analysis.

**Figure 9.** The dimensionless volume flow rate *Π<sup>δ</sup> V* (**a**) and dimensionless dissipation *Π<sup>δ</sup> Q* (**b**) as functions of the dimensionless pressure gradient *Π<sup>δ</sup> p*,*x* for a square-pitched screw with *t*/*D<sup>b</sup>* = **1.0** and a polymer melt with the power-law index *n* = 0.2.

#### *2.3. Data-Based Modeling*

#### 2.3.1. Symbolic Regression Analysis

Removing the need for numerical simulations required two analytical regressions that accurately predict the numerical solutions of our parametric design study. For this reason, we approximated the numerical solutions of our parametric design study analytically by using symbolic regression based on genetic programming implemented in the open-source software HeuristicLab (Hagenberg, Austria) [45]. This data-based modeling approach searches the space of mathematical expressions to find regressions that relate sets of input and output data. For detailed information, the reader is referred to [46,47].

To reduce the ranges of our target variables and the nonlinearities in our dataset, we took advantage of the following phenomenon: For large dimensionless pressure gradients, our numerical results showed that the influence of the screw-pitch ratio converges to zero, causing the pressure flow across the flight clearance to dominate the overall flow characteristics. Ignoring the effect of the flow in the *z*-direction, we used the following analytical relationships to approximate the dimensionless volume flow rate and the dissipation [26]:

$$\Pi\_{V,app}^{\delta} = 1 - \text{sign}\left(\Pi\_{p,\mathbf{x}}^{\delta}\right) \frac{\mathfrak{Z}^{\frac{1}{n}} n}{2n+1} \left|\Pi\_{p,\mathbf{x}}^{\delta}\right|^{\frac{1}{n}}.\tag{37}$$

$$
\Pi\_{Q,app}^{\delta} = \frac{\left(2n+1\right)^{n}}{n^{n}} \left| \Pi\_{V}^{\delta} - 1 \right|^{(1+n)}.\tag{38}
$$

Considering the flow of a power-law fluid, the first relationship describes the dimensionless volume flow rate as a function of the pressure gradient by a linear superposition of a one-dimensional drag and pressure flow. The second relationship, in contrast, approximates the dimensionless dissipation as a function of the volume flow rate by a one-dimensional pressure flow. These parameters were used to create new target variables by correcting the numerical results:

$$
\Delta \Pi\_V^\delta = \Pi\_V^\delta - \Pi\_{V,app}^\delta. \tag{39}
$$

$$
\Delta \Pi\_Q^\delta = \frac{\Pi\_Q^\delta}{\Pi\_{Q,app}^\delta + 1}.\tag{40}
$$

Figure 10 illustrates the characteristics of the modified target variables for a squarepitched screw and various power-law indices. Rather than showing a strong nonlinear increase for increased dimensionless pressure gradients, both parameters reach a plateau, thereby decreasing the value range of the target variables. Again, the curves are symmetrical about the point of pure drag flow (*Π<sup>δ</sup> <sup>p</sup>*,*<sup>x</sup>* = 0 or *Π<sup>δ</sup> <sup>V</sup>* = 1).

In the construction of the regression models, we randomly divided our dataset into three subsets: (i) a training set, (ii) a test set, and (iii) a validation set, including 4000, 2000, and 3231 design points, respectively. The first two subsets were employed to develop and optimize the symbolic regression solutions, while the third subset was used to validate the models against unseen numerical results not used in model development.

We applied the NSGA-II algorithm [48] to construct two regressions of the form:

$$
\Delta \Pi\_V^\delta = f \left( \Pi\_{p, \mathbf{x} \prime}^\delta \ n\_\prime \text{ t}/D\_\mathbf{b} \right); \tag{41}
$$

$$
\Delta \Pi\_Q^\delta = f \begin{pmatrix} \Pi\_{V'}^\delta & \mathfrak{n}\_\prime & \mathfrak{t}/\mathcal{D}\_\mathbb{b} \end{pmatrix}. \tag{42}
$$

This multi-objective non-dominant genetic algorithm simultaneously optimizes model quality and complexity. The latter was adjusted by restricting (a) the model size to a maximum tree length of 100, and (b) the function set, which defines the functions used to generate symbolic regression solutions, to (i) constant, (ii) state variable, (iii) addition, (iv) multiplication and division, and (v) cosine and sine, thus limiting the search space of the

regression analysis. Model optimization was driven by a constant optimization evaluator, which calculates Pearson *R* <sup>2</sup> of a symbolic regression solution and optimizes the constants used:

$$R^2 = 1 - \frac{\sum\_{i=1}^{n} (y\_i - \mathfrak{H}\_i)^2}{\sum\_{i=1}^{n} (y\_i - \overline{y})^2} \,\prime \tag{43}$$

where *y<sup>i</sup>* and *y*ˆ*<sup>i</sup>* are the numerical and approximated results, respectively, and *y* is the mean of the numerical solutions.

**Figure 10.** The corrected dimensionless volume flow rate **∆***Π<sup>δ</sup> V* (**a**) and corrected dimensionless dissipation **∆***Π<sup>δ</sup> Q* (**b**) as functions of the dimensionless volume flow rate *Π<sup>δ</sup> V* for a square-pitched screw with *t*/*D<sup>b</sup>* = **1.0** and a polymer melt with power-law index *n* = **0.2**.

For each model, we first performed 20 runs to generate a set of symbolic regression solutions, using the training and test data and then selected the most accurate approximation. To evaluate model quality, we carried out an error analysis for all subsets.

#### 2.3.2. Symbolic Regression Results

Our hybrid modeling approach provided two analytical regression for the corrected dimensionless volume flow rate ∆*Π<sup>δ</sup> V* and the corrected dimensionless dissipation ∆*Π<sup>δ</sup> Q* :

$$
\Delta\Pi\_V^\delta \left(\Pi\_{p,\mathbf{x}'}^\delta \ n\_\star \text{ t}/D\_b\right) = a\_{00} + \frac{A\_1 + A\_2}{A\_4 + A\_5 + A\_6} \text{;}\tag{44}
$$

$$
\Delta\Pi\_Q^\delta \left(\Pi\_{V'}^\delta \ n\_\prime \ t/D\_b\right) = b\_{00} + \frac{B\_1 + \frac{B\_2}{B\_3} + B\_4}{B\_5 + B\_6},\tag{45}
$$

where *A*1–*A*<sup>6</sup> and *B*1–*B*<sup>6</sup> are the subfunctions, which contain 25 (*a*00–*a*24) and 31 (*b*00–*b*30) coefficients, respectively. The subfunctions and their coefficients are given in the Appendix A. Note that the regressions predict the corrected target variables, while the actual parameters result from:

$$
\Pi\_V^\delta = \Delta \Pi\_V^\delta + \Pi\_{V,app}^\delta. \qquad \qquad \Pi\_Q^\delta = \Delta \Pi\_Q^\delta \left(\Pi\_{Q,app}^\delta + 1\right). \tag{46}
$$

Avoiding complex nested analytical functions, the models include only basic arithmetic operations and analytical functions. Due to their relatively simple mathematical structure, the approximations allow fast computation of the target variables without the need for further numerical simulations.

To evaluate the accuracy of the regression solutions, we carried out an error analysis for all subsets, including 9,231 design points in total. For this reason, we determined the mean absolute error (*AEmean*), the maximum absolute error (*AEmax*), and the coefficient of determination (*R* 2 ). The accuracy of the dissipation model was additionally investigated by comparing the mean relative error (*REmean*) and the maximum relative error (*REmax*):

$$AE\_{mean} = \frac{1}{N} \sum\_{i=1}^{n} |y\_i - \mathcal{Y}\_i|. \qquad \qquad AE\_{max} = \max\left(|y\_i - \mathcal{Y}\_i|\right). \tag{47}$$

$$RE\_{mean} = \frac{1}{N} \sum\_{i=1}^{n} \frac{|y\_i - \hat{y}\_i|}{y\_i}. \tag{28} \\ RE\_{max} = \max\left(\frac{|y\_i - \hat{y}\_i|}{y\_i}\right). \tag{48}$$

Overviews of the results of the error measures are given in Tables 5 and 6. Both models achieved a coefficient of determination *R* <sup>2</sup> > 0.999 for all subsets, which indicates excellent accuracy of the approximations. Since the results of the error measures of all subsets fall within a similar range, we conclude that overfitting was avoided. The low mean absolute and mean relative errors confirm the validity of the new models over the entire range of input parameters. Note that the design points of the validation set were not used for model development and therefore enabled an unbiased estimation of model quality. To demonstrate the performance of the regression models on the original numerical data, including 9231 design points, Table 7 shows the error measures of the final models in Equation (46). As the corrected target variables were constructed by using addition and multiplication, the models for the flow rate and dissipation exhibit the same absolute and relative errors, respectively. Comparisons of numerical and approximated results are illustrated in Figures 11 and 12 for various setups. The points indicate numerical results, and the continuous lines approximated solutions. Further, Figures 13 and 14 show a normalized representation of all design points in the form of scatter plots that compare numerical and approximated results for all subsets. Removing the need for numerical simulations, the new leakage-flow models can be used to include the effect of the flight clearance in the analysis of melt-conveying zones.

**Table 5.** The error measures of ∆*Π<sup>δ</sup> <sup>V</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> p*,*x* .


**Table 6.** The error measures of ∆*Π<sup>δ</sup> <sup>Q</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> V* .


**Table 7.** The error measures of *Π<sup>δ</sup> <sup>V</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> p*,*x* and *Π<sup>δ</sup> <sup>Q</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> p*,*x* 


.

**Figure 11.** A comparison of approximated results obtained from ∆*Π<sup>δ</sup> <sup>V</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> p*,*x* and numerical solutions for *t*/*D<sup>b</sup>* = 0.50 (**a**) and *t*/*D<sup>b</sup>* = 1.23 (**b**). The points indicate numerical results, and the continuous lines approximated solutions.

**Figure 12.** A comparison of approximated results obtained from ∆*Π<sup>δ</sup> <sup>Q</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> V* and numerical solutions for *t*/*D<sup>b</sup>* = 0.50 (**a**) and *t*/*D<sup>b</sup>* = 1.23 (**b**). The points indicate numerical results, and the continuous lines approximated solutions.

**Figure 13.** A scatter plot of ∆*Π<sup>δ</sup> <sup>V</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> p*,*x* : training and test set (**a**) and validation set (**b**). The dashed lines indicate an absolute error of 0.07.

**Figure 14.** A scatter plot of ∆*Π<sup>δ</sup> <sup>Q</sup>* = *f t*/*D<sup>b</sup>* , *n*, *Π<sup>δ</sup> V* : training and test set (**a**) and validation set (**b**). The dashed lines indicate a relative error of 10%.

#### **3. Conclusions**

We have proposed novel analytical regression models for predicting the volume flow rate and the viscous dissipation rate in the flight clearance of single-screw extruders. Unlike theories that correct the drag and pressure flows in the traditional pumping model, these approximations were designed to locally evaluate the characteristics of the leakage flow. Using a local reference system of the flight clearance allowed us to relax a variety of modeling assumptions and, therefore, to increase the accuracy of the models.

The approach presented here combines analytical, numerical, and data-based modeling techniques. The governing flow equations were rewritten in a dimensionless form by applying the theory of similarity. Three physically independent dimensionless input parameters of the flow were identified: (i) the screw-pitch ratio *t*/*D<sup>b</sup>* , (ii) the power-law index *n*, and (iii) a dimensionless pressure gradient across the flight clearance *Π<sup>δ</sup> p*,*x* . These dimensionless parameters were varied to create 9,231 physically independent setups, whose volume flow rates and dissipations were evaluated numerically with the shooting method. The numerical solutions of the design study were then approximated by means of symbolic regression based on genetic programming. The hybrid modeling approach yielded two

analytical relationships for volume flow rate and dissipation, the accuracy of which we have demonstrated in an error analysis, yielding a Pearson *R* <sup>2</sup> > 0.999 for both models. The first model showed a maximum absolute error of *AEmax* = 0.06735, while the second model produces a maximum relative error of *REmax* = 6.64%.

A novel feature of our theory is the flat-plate representation of the flight clearance. This local reference system gives rise to a two-dimensional flow in which the drag- and pressure flows in the cross- and down-channel directions are coupled due to the shearrate-dependent viscosity of the power-law fluid. To be able to solve the governing flow equations, we introduced a mathematical constraint by omitting the pressure gradient along the flight clearance (in the *z*-direction). From a physical viewpoint, this means that the corresponding flow component is governed exclusively by the rotation of the screw. This assumption was shown to be valid if the flight clearance is significantly smaller than the channel depth (*δ*/*h* 1).

To extend the applicability of the new leakage-flow models to a variety of screw designs and processing conditions in industrial use, our approach considered a significant range of dimensionless pressure gradients. Preliminary analyses indicated that these parameter ranges cover several orders of magnitude, depending on the depth of the leakage gap. Low values were observed for conventional screw designs with standard clearances (*δ* ≈ 0.001 *D<sup>b</sup>* ), while high values were detected for high-performance screws with pronounced undercut distances between main and secondary flight (e.g., wavedispersion screws). In addition, the calculations illustrated that the direction of leakage flow is governed mainly by the conveying characteristics of the extruder screw. For pressuregenerating melt-conveying zones, leakage flow reduces the net flow of the processing machine. For strongly overridden functional zones, in contrast, leakage flow increases the net flow.

Since the diverse characteristics of the leakage flow were considered in the construction of the design points, the numerical results of our parametric study extend over several orders of magnitude. To decrease the parameter ranges for the symbolic regression analyses, two modified target variables were constructed by assuming that, for large dimensionless pressure gradients, the flow can be approximated mathematically by linear superposition of one-dimensional drag and pressure flows. The corrected parameters were designed to reach a plateau for increased dimensionless pressure gradients, thereby decreasing the nonlinearities in the dataset.

The novel leakage-flow models will be implemented in our screw-simulation routine based on network theory (presented in [1–40]). The aim is to accurately describe the properties of the network elements positioned in the leakage gap to include the effect of the flight clearance in the prediction of pumping capability and power consumption. Using the new leakage-flow models in the network calculation, the modeling approach is no longer limited to any specific geometrical design and allows fast and accurate analysis of various types of extruder screws, including both conventional and high-performance screws.

**Author Contributions:** Conceptualization, C.M. and W.R.; Methodology, C.M. and W.R.; Software, C.M. and W.R.; Validation, C.M.; Formal Analysis, C.M.; Investigation, C.M. and W.R.; Resources, C.M., W.R., M.D., and G.S.; Data Curation, C.M.; Writing—Original Draft Preparation, C.M.; Writing— Review and Editing, W.R. and M.D.; Visualization, C.M. and W.R.; Supervision, W.R., M.D., and G.S.; Project Administration, C.M., W.R., M.D., G.S., and V.S.; Funding Acquisition, C.M, W.R., M.D., G.S, and V.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Austrian Science Fund (FWF): I 4872-N. This work was performed within the Competence Center CHASE GmbH, funded by the Austrian Research and Promotion Agency (grant number 868615). The authors acknowledge financial support by the COMET Centre CHASE (project No 868615), which is funded within the framework of COMET— Competence Centers for Excellent Technologies—by BMVIT, BMDW, and the Federal Provinces of Upper Austria and Vienna. The COMET program is run by the Austrian Research Promotion Agency (FFG).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Open Access Funding by the Austrian Science Fund (FWF).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Appendix A.1. Subfunctions of Equations (44) and (45)*

$$\left(\frac{t}{D\_b}\right)\_{eff} = \pi \tan\left(90 - \arctan\left(\frac{t}{D\_b\pi}\right)\right) \tag{A1}$$

$$A\_1 = a\_{01} \, II\_{p, \mathbf{x}}^{\delta} \left( \frac{t}{D\_b} \right)\_{eff} \tag{A2}$$

$$A\_2 = a\_{02} \left( a\_{03} + n \right) \left( a\_{04} + n \right) \Pi\_{p,x}^{\delta} \left( a\_{05} + \left( a\_{06} + n \right) \Pi\_{p,x}^{\delta} \left( + a\_{07} + a\_{08} \right) \left( \frac{t}{D\_b} \right)\_{eff} \right). \tag{A3}$$

$$A\_3 = n + \left(a\_{09} + a\_{10} \left. \Pi\_{p,x}^{\delta} \right. \right) \left(\frac{t}{D\_b}\right)\_{eff} \,. \tag{A4}$$

$$A\_4 = a\_{11} + a\_{12} \, \Pi\_{p, \text{x}}^{\delta \quad 2}.\tag{A5}$$

$$A\_5 = a\_{13} \left( a\_{14} + \Pi\_{p,x}^{\delta}{}^2 + a\_{15} \left( \frac{t}{D\_b} \right)\_{eff} \right) \left( a\_{16} \left( \Pi\_{p,x}^{\delta}{}^2 + \Pi\_{p,x}^{\delta}{}^4 + a\_{17} \left( \frac{t}{D\_b} \right)\_{eff} \right) . \tag{A6}$$

$$A\_{6} = a\_{18}n^{3}\left(a\_{19} + n + a\_{20}\left(\frac{t}{D\_{b}}\right)\_{eff}\right) + a\_{21}n^{2}\left(a\_{22} + n\right)\left(\Pi\_{p,x}^{\delta}{}^{2} + a\_{23}\left(\frac{t}{D\_{b}}\right)\_{eff}\right)\left(\Pi\_{p,x}^{\delta}{}^{2} + a\_{24}\left(\frac{t}{D\_{b}}\right)\_{eff}\right). \tag{A7}$$

$$B\_1 = b\_{01} + b\_{02}n + b\_{03} \left(\frac{t}{D\_b}\right)\_{eff} + b\_{04}n \left(\frac{t}{D\_b}\right)\_{eff} + n \left(\frac{t}{D\_b}\right)\_{eff} \left(b\_{05}n + b\_{06} \left(\frac{t}{D\_b}\right)\_{eff}\right) + b\_{07} \left(\frac{t}{D\_b}\right)\_{eff}^2. \tag{A8}$$

$$B\_2 = b\_{08} + b\_{09} \ II\_V^\delta + \left(\frac{t}{D\_b}\right)\_{eff} \left(b\_{10} + b\_{11} \ n^3\right) + \left(\frac{t}{D\_b}\right)\_{eff}^2 \left(b\_{12} + b\_{13} \ n^2\right). \tag{A9}$$

$$B\_3 = b\_{14} + b\_{15}n + b\_{16}\Pi\_V^\delta + b\_{17} \left(\Pi\_V^\delta\right)^2 + \frac{1}{b\_{18} + b\_{19}\left(\Pi\_V^\delta\right)^2} + b\_{20}\left(\frac{t}{D\_b}\right)\_{eff}.\tag{A10}$$

$$B\_4 = b\_{21} \sin \left( b\_{22} \, \Pi\_V^\delta \right). \tag{A11}$$

$$B\_5 = b\_{23} + b\_{24} \, \Pi\_V^\delta + b\_{25} \left(\Pi\_V^\delta\right)^2 + b\_{26} \left(\frac{t}{D\_b}\right)\_{eff}.\tag{A12}$$

$$B\_6 = n \left( b\_{27} + \frac{b\_{28}}{b\_{29} - \left( I I\_V^\delta \right)^2} + b\_{30} \left( \frac{t}{D\_b} \right)\_{eff} \right). \tag{A13}$$

#### *A.2. Model Coefficients*

#### **Table A1.** Model coefficients.



**Table A2.** Nomenclature.

#### **References**

2. White, J.L.; Potente, H. *Screw Extrusion*; Hanser Publishers: Munich, Germany, 2001; ISBN 9783446196247.

<sup>1.</sup> Tadmor, Z.; Klein, I. *Engineering Principles of Plasticating Extrusion*; Van Nostrand Reinhold, Co.: New York, NY, USA; ISBN 9780442156350.

