Applying network theory to the modeling of multilayer flows in slot dies: a use case for symbolic regression-based co-extrusion prediction models

1 Introduction

Extruded films, sheets, semi-finished sheets for thermoforming, and coatings for substrate webs are essential not only in packaging applications, but also in various other industry sectors, such as agriculture, automotive, healthcare, and consumer goods. These products – characterized by their large width-to-height ratios – are commonly manufactured by extrusion through slot dies (i.e., dies with slit-shaped exit cross-sections) that transform initially round melt strands into flat, rectangular extrudates. The distribution of melt across the width is generally achieved by a so-called manifold section, which is followed by a flow-resistance region that is commonly referred to as die land (Hopmann and Michaeli 2016). In all types of plastics processing technology, including the extrusion of flat films and sheets, both product quality and efficiency of the related manufacturing process are determined mainly by the material properties, the processing conditions, and the geometry of the processing equipment. Consequently, in melt extrusion (i) choosing polymers with suitable rheological properties, (ii) optimizing the design of the flow domain of the shaping die (i.e., the slot die), and (iii) choosing proper processing parameters, such as throughputs and temperatures, are key to a successful process, especially when different polymer types are extruded simultaneously to create multilayer structures with properties that are superior to those of the individual components. In these co-extrusion techniques, so-called feedblock systems are often used to layer individual melt streams on top of each other and guide the stratified flow into a conventional slot die for the final shaping operation. Hatzikiriakos and Migler (2005) summarized common challenges in flat-film and sheet extrusion: (i) melt-fracture effects (such as shark skin) at high production rates, (ii) neck-in effects, and – particularly in co-extrusion – (iii) layer encapsulation and (iv) interfacial flow instabilities. Furthermore, (v) die deflection under operating pressure (clam shelling) can lead to a non-uniform thickness distribution across the width (Deshmukh et al. 2020; Helmy 1988; Pfeiffer and Dietz 1997).

Optimizing die and process designs usually requires solving a mathematically inverse problem (Elgeti et al. 2012; Hao et al. 2023; Yilmaz and Kirkkopru 2020). Hence, full factorial or optimization algorithm-supported parameter studies are often necessary to find optimal designs and conditions. In this context, it is advantageous to apply modeling and simulation techniques that are (i) able to accurately predict flow parameters (e.g., pressure-throughput behavior and shear stresses) that determine quality and efficiency, (ii) fast, (iii) mathematically stable, and (iv) easy to parametrize (generalize). The semi-numerical network theory approach (also called electrical network theory) has evolved in computational flow modeling as an alternative to (i) time- and resource intensive numerical methods (e.g., computer fluid dynamics (CFD) simulations), (ii) analytical methods that oversimplify, are inaccurate or are not applicable, and to (iii) purely data-driven methods that do not take into account the actual physical phenomena of a process. In network theory the analogy between electric circuit and fluid dynamics is used to solve such flow problems fast and efficiently while fulfilling the aforementioned requirements. Consequently, network theory has been used to model the polymer melt flow and temperature behavior in single-screw extruders (e.g., Luger and Miethlinger 2024; Marschik et al. 2020; Roland et al. 2020; Zacher et al. 2023) and extrusion dies (Booy 1982; Köpplmayr and Miethlinger 2013; Razeghiyadaki et al. 2022).

Here, we present an efficient and accurate network theory-based approach to modeling the pressure-throughput behavior of non-Newtonian two-layer co-extrusion flows through slot dies that have a rectangular manifold channel. The remainder of this article is structured as follows: First, state-of-the-art modeling techniques for slot dies and their limitations are briefly discussed, followed by a summary of the hybrid modeling approach used to develop the symbolic regression models for predicting one-dimensional and isothermal two-layer co-extrusion flows of non-Newtonian fluids through simple rectangular ducts and how they are integrated into network theory. Finally, the results from network theory are presented and compared to the results of CFD simulations with a commercial simulation software package.

2 State-of-the-art modeling of slot dies

The main purpose of a slot die is to shape an incoming melt stream according to the requirements of a downstream processing step (e.g., stretching and cooling on a chill roll or shaping in a polishing stack). Ideally, the flow in the slot die will (i) be distributed across the width homogeneously in order to produce a film of constant thickness, (ii) have a highly uniform history of temperature, residence time, and level of deformation along all flow paths, and (iii) meet requirements (i) and (ii) for a wide range of processing conditions and material properties.

A large variety of slot die designs have been established, which can be classified mainly by the cross-sectional shape of the manifold channel (e.g., rectangular, circular, and tear-drop shaped) and the length distribution of the subsequent die land across the die width (constant, linear, and non-linear). The simplest approach commonly used to create the initial die design is analytical flow analysis (Hopmann and Michaeli 2016). Here, a linear reduction of the flow rate across the width of the manifold channel and equal flow resistances (i.e., pressure drops) along all flow paths are considered. Assuming Newtonian flow behavior in a representative viscosity approach and equal shear rates both in the manifold and the land yields a design that is independent from operating conditions. For the mono-extrusion flow through a slot die with a slit-shaped manifold with constant manifold width w, the manifold height h(x) is:

and the die-land length y(x)

where x is the width coordinate with its origin at the lateral edge of the die and δ is the die land height. To obtain more accurate results and avoid computationally expensive and time-consuming numerical methods (e.g., CFD simulations), an electrical network theory approach can be employed. Network theory divides an initially complex flow domain into small segments of constant cross-sectional shape (i.e., sub-elements) whose pressure-throughput behavior (i.e., the characteristic die line) can be modeled analytically. This relationship is generally defined by the die conductance k, which depends on geometry parameters, processing conditions, and the rheological behavior of the fluid in a particular sub-element.

To capture the rheological behavior of the fluid, either a Newtonian model with a constant viscosity η or a power-law formulation with consistency K and the power-law index n can be implemented:

In both cases, use of a more complex rheological model from which a constant representative viscosity or a local power law can be derived is recommended. This increases modeling accuracy while simultaneously low mathematical complexity. For example, the five-parameter Carreau-Yasuda model (Yasuda 1979) captures the overall viscosity behavior of polymer melts across a wide range of shear rates:

The viscosities at zero shear rate and infinite shear rate are denoted by η₀ and η_∞, respectively; n_cy is the Carreau-Yasuda power-law index, λ the relaxation time, and γ ˙ the shear rate. The breadth of the transition region is defined by parameter a, and the temperature dependency of the viscosity is considered in the temperature shift factor a _T . The temperature dependence can be modeled by a variety of approaches, for instance, an approximate Arrhenius relationship:

where α is the temperature sensitivity of the viscosity and T₀ is the reference temperature. As shown schematically in Figure 1, Eq. (4) can be evaluated at a local effective shear rate γ ˙ r e p (which is representative of a given sub-element) to obtain a representative viscosity η_rep. Alternatively, the parameters K and n of a local power law, which can be interpreted as a local tangent to the Carreau-Yasuda model, can be determined by (as proposed for die flows by Marschik et al. [2019]; Marschik et al. [2020]):

Figure 1:

Determining a representative Newtonian viscosity (blue) and the parameters of a local power law (red) based on a representative (local effective) shear rate and the Carreau-Yasuda viscosity model (black).

The local effective (representative) shear rate for a simple mono-extrusion flow in a rectangular slit is a function of the geometry of the local element, the processing conditions (flow rate V ˙ ), and the power-law index n:

The pressure-throughput behavior of the sub-elements, which are now well-defined by aforementioned rheological behavior and geometry parameters (provided by the user), is given by:

Note that Eq. (9) can be simplified for Newtonian flow behavior by assuming:

The elements can then be interconnected to create a flow-resistance network. In analogy to electrical circuits, pressure differences can be expressed as voltages, while volume flow rates can be considered as currents. Applying Kirchhoff’s laws yields a system of equations that can be solved numerically by means of a mesh or a nodal analysis.

Note that, up to this point, both the analytical and network theory-based approach are strictly limited to mono-extrusion slot-die flows. When multilayer flows (i.e., co-extrusion processes) are to be considered, explicit expressions for the die conductance in a slit K′ are not available, iterative or numerical approaches (e.g., based on finite-element methods) using bespoke solvers or commercial CFD software packages are required. Hence, this work introduces symbolic regression models predicting co-extrusion flows through simple rectangular ducts into the electrical network theory approach, as discussed in detail the following two sections.

3 Hybrid modeling approach for co-extrusion flows

In previous studies, we derived (Hammer et al. 2021a, 2021b) and validated (Hammer et al. 2022) symbolic regression models that predict the one-dimensional and isothermal two-layer co-extrusion flow of non-Newtonian fluids through rectangular ducts. The hybrid modeling approach – using analytical, numerical, and data-based modeling steps – is briefly summarized below: assuming an isothermal, laminar, fully developed, stationary, one-dimensional, and stratified flow of two incompressible and wall-adhering power-law fluids through a rectangular duct of infinite width (the retarding effect of lateral walls is ignored), the momentum equation in the down-channel direction can be simplified to:

where ∂p/∂z is the down-channel pressure gradient, y the channel-height coordinate, v _z the down-channel velocity, and K _i and n _i the power-law parameters of fluids A and B (indication by index i). Solving Eq. (11) for both fluids A and B requires coupling of the results via the interfacial conditions that specify equal shear stresses and flow velocities at the interfacial position. This flow problem is then transformed into dimensionless representation by properly scaling the physical parameters to obtain two dimensionless fluid layers of unit height. As a result, rather than by seven parameters in dimensional representation, this co-extrusion flow problem can be fully described by four dimensionless input parameters: (i) the power-law index of the bottom-layer fluid A n _A , (ii) the power-law index of the top-layer fluid B n _B , (iii) the dimensionless position of the interface κ or the dimensionless volume flow rate of fluid A Π V A ,

and (iv) the ratio of dimensionless pressure gradients χ,

with h _A denoting the height of the bottom layer fluid A, h the overall channel height, V ˙ A and V ˙ B the individual volume flow rates of fluids A and B, respectively, and v_ref the mean flow velocity within the channel. Various dimensionless target quantities were defined (Hammer et al. 2021a, 2021b), where (i) the dimensionless pressure gradient Π p , z A and (ii) the dimensionless average shear rates for layer A γ ‾ ˙ A * and layer B γ ‾ ˙ B * , respectively, are relevant to this work:

The volumetric flow rate per unit of width can be determined by integration of the velocity profile over the channel height. For solving Eq. (11) in dimensionless representation (i.e., computation of the dimensionless velocity profile and the target quantities derived thereof), the numerical shooting method is employed. This method is perfectly suited to this kind of mathematical problem and allows a full-factorial design study that encompasses more than 44,000 design points to be evaluated within less than an hour on a desktop computer. This study was performed by varying the four independent influencing parameters over wide ranges to cover almost all combinations of material parameters, geometry, and process conditions that are practically relevant in co-extrusion.

The results of the numerical design study were then approximated by means of symbolic regression modeling based on genetic programming using the open-source environment HeuristicLab (Wagner et al. 2014). Using only very simple mathematical operators, this data-based modeling approach yields simple algebraic relationships between each of the target quantities and the four dimensionless input parameters:

Note that the detailed model that predicts the dimensionless pressure gradient Π p , z A has already been published (Hammer et al. 2021b), whereas the models for the dimensionless average shear rates γ ‾ ˙ A * and γ ‾ ˙ B * were derived specially for this work and can be found in the Supplementary Material. The main advantages of these models are their execution speed (i.e., real-time capability), prediction accuracy (since the physics of the underlying process is considered), numerical stability, and generalizability to a wide range of values of the input parameters. This hybrid modeling approach was discussed in more detail by Hammer et al. (2021b).

4 Network analysis based on symbolic regression co-extrusion models

Integration of the symbolic regression models into network theory to predict the two-layer co-extrusion flow behavior through slot dies using a numerical procedure to solve the network was realized in an in-house simulation routine implemented in Python. The overall workflow of the program is summarized in the flow chart in Figure 2.

Figure 2:

Flow chart of the simulation routine.

First, the input parameters for the simulation are defined, which include:

1. material properties structured in a database:

   1. material name,

   2. the parameters of the Carreau-Yasuda model,

   3. the reference temperature and the temperature sensitivity of the approximate Arrhenius relationship to evaluate the temperature shift factor a _T , and

   4. either a constant melt density or the parameters of the Tait equation;

2. process parameters:

   1. the individual mass throughputs of bottom-layer fluid A and top-layer fluid B,

   2. assignment of materials to fluids A and B (layer arrangement), and

   3. their respective melt temperatures;

3. geometry parameters:

   1. the half-width of the slot die W ,

   2. the manifold width w ,

   3. a mathematical function for the relationship between manifold height and the width-coordinate h(x),

   4. the die-land gap height δ, and

   5. a mathematical function for the relationship between die-land length and the width coordinate y(x);

4. numerical settings:

   1. the number of elements per half-width of the die,

   2. convergence criteria for the iterative procedure of the local power-law evaluation and the network analysis itself.

According to the input data, the flow domain is discretized to create the flow-resistance network under consideration of the symmetry condition at the center position, as illustrated in Figure 3 for a small number of elements. Next, the local properties for each sub-element are determined, which includes four major steps: (i) determination of the local power-law parameters according to Eqs. (6) and (7) for both fluids A and B, (ii) computation of the additional dimensionless parameters χ and Π V A using Eqs. (12) and (13), respectively, (iii) evaluation of the symbolic regression models for the dimensionless average shear rates γ ‾ ˙ A * and γ ‾ ˙ B * , and (iv) back-transformation of the dimensionless average shear rates into dimensional form (Eqs. (16) and (17)) to be used as representative shear rate. Steps (i) to (iv) are carried out iteratively until a converged solution for the power-law parameters is obtained. Note that for the first run of this “inner iteration”, the local effective shear rates for both fluid layers are initialized for each sub-element by the following simple definition:

Figure 3:

Schematic representation of the transformation of the slot die geometry into a flow resistance network.

The iterative procedure (“outer iteration loop”) for solving the network itself first requires the evaluation of the die conductances. Given that electrical network analysis of hydraulic systems requires simple analytical relationships between pressure drop and volume flow rate (i.e., the characteristic die line) for each of the sub-elements, the die conductance K', which depends on geometry parameters and rheological properties of the material, must be calculated. These simple analytical relationships are available for mono-extrusion flows, but not for co-extrusion flows. The hybrid modeling approach yields a simple but accurate prediction model for the dimensionless pressure gradient in down-channel direction Π p , z A with respect to the bottom fluid A. This dimensionless pressure gradient can be transformed back into dimensional form, that is ∂ p ∂ z , by Eq. (22). Due to the similarity of this equation to the characteristic die line of the mono-extrusion flow of a power-law fluid through slits:

we can derive an analytical expression for the die conductance, denoted as the equivalent die conductance for two-layer co-extrusion flows K'':

Consequently, the four dimensionless input parameters are again used to evaluate the dimensionless down-channel pressure gradient and subsequently calculate K''. Each of the sub-elements in Figure 3 actually consists of a parallel connection of a resistance and a source, as presented by Marschik et al. (2019). The first represents the pressure flow m ˙ p , i and the latter the theoretical drag flow m ˙ d , i of the local die segment, which is necessary to linearize the characteristic die line of non-Newtonian flows around the operating point. The total mass flow rate of a sub-element m ˙ i is thus modeled by:

where k _i is the linearized die conductance (i.e., the local slope of linearization) and p_in,i and p_out,i the pressure values at the neighboring nodal points. Note that the theoretical drag flow is not of physical relevance and represents the intercept of the linearized die line with y-axis (mass flow). Applying Kirchhoff’s current law, stating that the sum of all incoming and outgoing currents (mass flows) is equal to zero for each nodal point:

the network equations can be formulated according to the schematic representation in Figure 4. For an arbitrary node i one obtains:

where the subscripts m and l indicate the elements from manifold and land, respectively. The boundary nodes, i.e., the first node (i = 0) at width position x = 0 and last node (i = n) at the center position of the die x = W are considered as follows:

Figure 4:

Network structure for nodal points zero, i, and n.

Summarizing these equations for all nodes, a linear set of equations can be formulated:

In this framework, m ˙ includes the boundary conditions for the mass flow rate at the inlet and the outlet pressure, m ˙ d is the vector of drag flow rates, k is the matrix of linearized conductances, and p is the pressure field. Rearranging Eq. (29) by calculating the inverse of the conductance matrix, the network can be solved for the pressure field:

At the end of each iteration, the updated pressure field is used to evaluate the mass flow rates for each of the sub-elements. This outer iterative procedure is considered converged as soon as the user-defined convergence level for the difference of the pressure at the die inlet is reached.

The final post-processing step is to determine the overall pressure drop and the distribution of the mass-flow rate across the width of the slot die.

5 Simulation setup

A simple coat-hanger type laboratory-scale slot die was used for modeling. The manifold width w was considered constant, while both the manifold height h and the length of the die land y decreased according to a square-root function from the center towards the die edges, as summarized in Table 1. The half-width of the die was discretized into 2000 elements.

Table 1:

Geometry parameters of the coat-hanger type slot die modeled.

Geometry parameter	Value	Unit	Schematic representation
Manifold width w	15	mm
Manifold height h(x)	δ · x / w	mm
Die half width W	60	mm
Land length y(x)	2 x · w	mm
Land gap height δ	1	mm

A two-layer co-extrusion flow of commercially available extrusion grades was modeled in the network analysis. Acrylonitrile butadiene styrene copolymer (ABS) and polymethyl methacrylate (PMMA) were defined by a constant melt density. The shear-rate and temperature dependencies of the viscosities were described by a Carreau-Yasuda and an approximate Arrhenius model, respectively (see Table 2 and Figure 5). The simulations were conducted under isothermal conditions with melt temperatures of 240 °C. The total mass-flow rate was equally distributed between both materials and varied across a wide range to obtain different levels of the representative shear rates for the local power-law formulation. The numerical values of the material and process parameters are summarized in Table 2.

Figure 5:

Carreau-Yasuda models at 240 °C for ABS (solid line) and PMMA (dashed line) and values measured by a high-pressure capillary (high shear rates) and a plate-plate rheometer (shear rates smaller than 1 s⁻¹), respectively.

To validate the results of applying network theory, CFD simulations of the full three-dimensional flow domain were performed using the commercial software environment Ansys Fluent 2023 (ANSYS Inc., Canonsburg, PA, USA), which is based on the finite-volume method (FVM). The mesh size was approximately 4.5 million cells including an inflation layer with five cells towards the die walls and across the interface. To capture the co-extrusion flow behavior accurately, a multiphase configuration based on the volume-of-fluid method was used. For both the network-theory approach and the CFD simulations, a symmetry condition at the center position of the die was considered.

6 Results and discussion

The two most important parameters characterizing the flow in slot dies for films and sheets were evaluated: (i) the overall pressure consumption and (ii) the mass-flow distribution across the outlet of the die. Figure 6 shows a comparison of the pressure drops predicted by the network theory-based approach and the CFD simulations. The two methods are in good agreement across all mass-flow rates under investigation, as indicated by a mean relative error of approximately 4.4 % and a maximum relative error of 6.5 %. To assess the die’s ability to balance the material flux, the mass throughput per unit of width was evaluated as a function of the die half width W. To capture the full range of simulations S1 to S5, a logarithmic scale was chosen for the throughput (y-axis). For the CFD simulations, the throughput per unit of width was evaluated at distinct positions: 0.05 mm, 5 mm, 10 mm, 20 mm, 30 mm, 40 mm, 50 mm, and 60 mm. Again, network theory was capable of accurately predicting the general trend of a throughput distribution that is relatively uniform across the majority of the width and shows a strong gradient tending to zero towards the outer die wall (position x = 0). Disregarding the prediction at position 0.05 mm – where the greater deviation can be explained by the strong gradient and thus mesh dependency – a mean relative error between network and CFD simulations of approximately 3.3 % was obtained.

Figure 6:

Comparison between the results obtained from network theory and from CFD simulations: (a) pressure drop and (b) mass throughput per unit of width across the die width.

Minor disagreements between the results of the two approaches are most likely due to the following differences in modeling assumptions:

1. A local power-law formulation was used in network theory, while the flow behavior in the CFD simulations was captured by a more comprehensive Carreau-Yasuda model.

2. The symbolic regression models implemented in the network do not take into account the retarding effect of the lateral walls and also exhibit a slight deviation from underlying simulation data. This particularly influences the prediction accuracy of the regression models in regions of the manifold with height-to-width ratios greater than 0.1.

3. Network theory neither considers any curvature of the manifold nor transverse flows.

4. Network theory assumes a constant layer distribution in the die.

In summary, our approach provided highly accurate results within seconds on a notebook PC (HP EliteBook 845, with 4.0 GHz, and 32 GB RAM) and thus clearly outperformed the CFD simulation (duration of approximately two to three hours).

7 Conclusions & outlook

In this work, we have demonstrated the potential of applying symbolic regression models predicting the co-extrusion flow in simple ducts within a more global modeling framework. For this purpose, the pressure-throughput behavior of the two-layer co-extrusion flow through slot dies with a rectangular manifold channel was considered using a semi-numerical network theory-based approach. Symbolic regression models were employed to derive an equivalent die conductance factor and evaluate the average shear rate per layer (in order to implement a local power law) for each element of the network. A comparison to CFD simulation confirmed that our method not only offers high prediction accuracy, but is also computationally efficient, and thus constitutes a promising approach to finding optimal die geometries. Future activities will address further development of this network-theory approach with regard to: (i) the evaluation of additional flow parameters in the slot die (such as layer distribution, interfacial shear stress, residence time, and viscous dissipation), (ii) consideration of non-isothermal flow behavior, (iii) a higher degree of freedom in die geometry, and (iv) modeling of multilayer flows with an arbitrary number of layers.