Boundary-Conforming Finite Element Methods for Twin-Screw Extruders: Unsteady - Temperature-Dependent - Non-Newtonian Simulations

TL;DR

This paper introduces a boundary-conforming space-time finite element method with a novel mesh update technique for efficient, accurate simulations of unsteady, temperature-dependent non-Newtonian flows in twin-screw extruders, validated by experiments.

Contribution

The paper presents the SRMUM mesh update method enabling efficient, re-mesh-free simulations of complex twin-screw extruder flows with temperature dependence.

Findings

Mesh convergence demonstrated in 2D and 3D cases.

Agreement with experimental results confirmed.

Method successfully applied to complex industrial scenarios.

Abstract

We present a boundary-conforming space-time finite element method to compute the flow inside co-rotating, self-wiping twin-screw extruders. The mesh update is carried out using the newly developed Snapping Reference Mesh Update Method (SRMUM). It allows to compute time-dependent flow solutions inside twin-screw extruders equipped with conveying screw elements without any need for re-meshing and projections of solutions - making it a very efficient method. We provide cases for Newtonian and non-Newtonian fluids in 2D and 3D, that show mesh convergence of the solution as well as agreement to experimental results. Furthermore, a complex, unsteady and temperature-dependent 3D test case with multiple screw elements illustrates the potential of the method also for industrial applications.

Tables

Table 1. Table 1: Screw geometry parameters

Screw radius $R_s$	15.275 $mm$
Center line distance $C_l$	26.2 $mm$
Screw-screw clearance ${\delta }_s$	0.2 $mm$
Screw-barrel clearance ${\delta }_b$	0.15 $mm$

Table 2. Table 4: Geometry parameters for 3D forward-conveying screw element.

Screw radius $R_s$	14.75 $mm$
Center line distance $C_l$	26.5 $mm$
Screw-screw clearance ${\delta }_s$	0.25 $mm$
Screw-barrel clearance ${\delta }_b$	0.25 $mm$
Pitch length	28 $mm$

Table 3. Table 6: Geometry parameters for 3D forward-conveying screw element used in convergence study for shear-thinning fluid.

Screw radius $R_s$	16.00 $mm$
Center line distance $C_l$	25.5 $mm$
Screw-screw clearance ${\delta }_s$	0.25 $mm$
Screw-barrel clearance ${\delta }_b$	0.25 $mm$
Pitch length	28 $mm$

Table 4. Table 8: Geometry parameters of a 3D conveying screw element for temperature-dependent flow.

Screw radius $R_s$	156.00 $mm$
Center line distance $C_l$	262.0 $mm$
Screw-screw clearance ${\delta }_s$	4.0 $mm$
Screw-barrel clearance ${\delta }_b$	4.0 $mm$
Pitch length $p_l$	280 $mm$

Table 5. Table 11: Mesh discretization for 3D unsteady temperature-dependent flow.

mesh	$n_s$	$n_r$	$n_a$	$\#$ elements
1	180	10	450	1 620 000
2	360	20	900	12 960 000

Equations

$${\mathbf{x}}_i={\mathbf{x}}_{screw}(I{D}_i)+d_{rel}(I{D}_i)\ast \left[{\mathbf{x}}_{barrel}(I{D}_i)-{\mathbf{x}}_{screw}(I{D}_i)\right].$$

$$y_{rel}=C{P}_{rel}=ID/I{D}_{CP}.$$

$$y_{midLine}=C{P}_{rel}\ast R_{barrel}\ast sin(C{P}_{rel}\ast {\theta }_{CP})+(1-C{P}_{rel})\ast avg(y_{screw}),$$

$$\nabla \cdot \mathbf{u}=0\text{mbox}on{\Omega }_t,\forall t\in (0,T),$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)-\nabla \cdot \mathbf{\sigma }=0\text{mbox}on{\Omega }_t,\forall t\in (0,T),$$

$$\rho c_p\left(\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T\right)-\kappa \mathbf{\Delta }T-\underset{\text{viscous dissipation}\varphi }{\underbrace{2\eta \nabla \mathbf{u}:\mathbf{\epsilon }\left(\mathbf{u}\right)}}=0\text{mbox}on{\Omega }_t,\forall t\in (0,T),$$

$$\mathbf{\sigma }(\mathbf{u},p)=-p\mathbf{I}+2\eta \left(\dot{\gamma },T\right)\mathbf{\epsilon }(\mathbf{u}),$$

$$\mathbf{\epsilon }(\mathbf{u})=\frac{1}{2}\left(\nabla \mathbf{u}+\left({\nabla \mathbf{u}}^T\right)\right).$$

$$\mathbf{u}={\mathbf{g}}^f\text{mbox}on{\left({\Gamma }_t\right)}_g^f,$$

$$\mathbf{n}\cdot \mathbf{\sigma }={\mathbf{h}}^f\text{mbox}on{\left({\Gamma }_t\right)}_h^f,$$

$$T=g^T\text{mbox}on{\left({\Gamma }_t\right)}_g^T,$$

$$\mathbf{n}\cdot \kappa \nabla T=h^T\text{mbox}on{\left({\Gamma }_t\right)}_h^T.$$

$$\dot{\gamma }=\sqrt{2\mathbf{\epsilon }\left(\mathbf{u}\right):\mathbf{\epsilon }\left(\mathbf{u}\right)}.$$

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

$$\eta \left(\dot{\gamma },T\right)=\frac{{\eta }_0\left(T\right)}{1+{\left(\frac{{\eta }_0\left(T\right)\dot{\gamma }}{{\tau }^{\ast }}\right)}^{(1-n)}},$$

$${\eta }_0(T)=D_{1}exp\left(-\frac{A_1\left(T-T_{ref}\right)}{A_2+\left(T-T_{ref}\right)}\right),$$

$$({\mathcal{S}}_{\mathbf{u}}^h{)}_n=\{{\mathbf{u}}^h\in [H^{1h}(Q_n){]}^{n_{sd}}|{\mathbf{u}}^h\doteq {\mathbf{g}}^{f,h}\text{on}(P_n{)}_{\mathbf{g}}\},$$

$$({\mathcal{V}}_{\mathbf{u}}^h{)}_n=\{{\mathbf{w}}^h\in [H^{1h}(Q_n){]}^{n_{sd}}|{\mathbf{w}}^h\doteq 0\text{on}(P_n{)}_{\mathbf{g}}\},$$

$$({\mathcal{S}}_p^h{)}_n=({\mathcal{V}}_p^h{)}_n=\{p^h\in H^{1h}(Q_n)\},$$

$$({\mathcal{S}}_T^h{)}_n=\{T^h\in H^{1h}(Q_n)|T^h\doteq g^{T,h}\text{on}(P_n{)}_g\},$$

$$({\mathcal{V}}_T^h{)}_n=\{v^h\in H^{1h}(Q_n)|v^h\doteq 0\text{on}(P_n{)}_g\}.$$

$$\begin{array}{r}{\int }_{Q_n}{\mathbf{w}}^h\cdot \rho \left(\frac{\partial {\mathbf{u}}^h}{\partial t}+{\mathbf{u}}^h\cdot {\nabla \mathbf{u}}^h\right)dQ+{\int }_{Q_n}\mathbf{\epsilon }({\mathbf{w}}^h):{\mathbf{\sigma }}^h(p^h,{\mathbf{u}}^h)dQ\\ +{\int }_{Q_n}{q}^h\nabla \cdot {\mathbf{u}}^{h}dQ+{\int }_{{\Omega }_n}({\mathbf{w}}^h{)}_n^{+}\cdot \rho \left(({\mathbf{u}}^h{)}_n^{+}-({\mathbf{u}}^h{)}_n^{-}\right)d\Omega \\ +\sum _{e=1}^{({n_{el})}_n}{\int }_{Q_n^e}{\tau }_{\text{mbox}MOM}\frac{1}{\rho }\left[\rho \left(\frac{\partial {\mathbf{w}}^h}{\partial t}+{\mathbf{u}}^h\cdot {\nabla \mathbf{w}}^h\right)\right]\\ \cdot \left[\rho \left(\frac{\partial {\mathbf{u}}^h}{\partial t}+{\mathbf{u}}^h\cdot {\nabla \mathbf{u}}^h\right)-\nabla \cdot {\mathbf{\sigma }}^h(p^h,{\mathbf{u}}^h)\right]dQ\\ +\sum _{e=1}^{({n_{el})}_n}{\int }_{Q_n^e}{\tau }_{\text{mbox}CONT}\nabla \cdot {\mathbf{w}}^h\rho \nabla \cdot {\mathbf{u}}^{h}dQ={\int }_{P_n}{\mathbf{w}}^h\cdot {\mathbf{h}}^{f,h}dP.\end{array}$$

$$\begin{array}{r}{\int }_{Q_n}{v}^h\cdot \rho c_p\left(\frac{\partial T^h}{\partial t}+{\mathbf{u}}^h\cdot \nabla T^h\right)dQ+{\int }_{Q_n}\nabla v^h\cdot \kappa \nabla T^{h}dQ\\ -{\int }_{Q_n}{v}^h\varphi dQ+{\int }_{{\Omega }_n}(v^h{)}_n^{+}\rho c_p\left((T^h{)}_n^{+}-(T^h{)}_n^{-}\right)d\Omega \\ +\sum _{e=1}^{({n_{el})}_n}{\int }_{Q_n^e}{\tau }_{\text{mbox}TEMP}\frac{1}{\rho c_p}\left[\rho c_p\left(\frac{\partial v^h}{\partial t}+{\mathbf{u}}^h\cdot \nabla v^h\right)\right]\\ \cdot \left[\rho c_p\left(\frac{\partial T^h}{\partial t}+{\mathbf{u}}^h\cdot \nabla T^h\right)-\nabla \cdot \kappa \nabla T^h-\varphi \right]dQ\\ ={\int }_{P_n}{v}^h{h}^{T,h}dP.\end{array}$$

$$\begin{array}{r}{\left({\mathbf{u}}^h\right)}_n^{\pm }=\underset{\zeta \to 0}{\lim }{\mathbf{u}}^h\left(t_n\pm \zeta \right)\\ {\int }_{Q_n}...dQ={\int }_{I_n}{\int }_{{\Omega }_n}...d\Omega dt\\ {\int }_{P_n}...dP={\int }_{I_n}{\int }_{{\Gamma }_n}...d\Gamma dt\end{array}$$

Full text

††journal: Computers $\&$ Fluid

1 Introduction

Co-rotating twin-screw extruders are very important processing devices within the plastic-producing industry. A single twin-screw extrusion process allows to carry out multiple operations at the same time - melting, compounding, blending, pressurization, shaping. Furthermore, the modular structure of the extruder enables to design configurations that are tailored to specific processes. This meets the need of the industry for more sophisticated and specialized polymers. Typical screw elements are conveying, kneading and mixing elements. Twin-screw extruders are used in particular for compounding processes since they provide short residence times and extensive mixing. In order to improve new screw configurations and by that new processes, predictions of residence time, mixing behavior as well as temperature distributions are necessary. This requires detailed information about the flow inside the extruder. Obtaining this information through experiments is a time-consuming and difficult task, if not impossible. This is due to the complex moving domain, small gap sizes and high pressures. Furthermore, experiments are not able to provide multiple classes of information - like viscosity distributions, temperature and pressure - at a time. This makes the computation of the flow inside twin-screw extruders using Computational Fluid Dynamics (CFD) extremely promising and appealing. By that, different screw desings, materials and flow conditions can be tested at the same time by running different simulations. Furthermore, the flow results can be post-processed in many different ways, allowing to extract various quantities of interest (Mierka et al., 2014).

There exist a broad range of models trying to simulate the flow inside twin-screw extruders. First results have been obtained in (Chen and White, 1991) using a 1D model. It allows to predict the pressure and temperature build-up along the axial length of the extruder. 2D results using a fictious domain method with mesh refinement were presented in (Bertrand et al., 2003). However, a full 3D simulation is necessary to obtain detailed information about the flow field. The major challenge for simulations in twin-screw extruders is the complex, constantly changing domain with small gap sizes, making meshing extremely challenging. Several numerical methods have been used to tackle this issue:

A mesh superposition method that avoids remeshing has been incoporated into Polyflow (Zhang et al., 2009). Individual meshes are created for each screw and the related flow domain and then coupled by interpolation. In (Zhu et al., 2013), this method was used to compute the flow field inside tri-screw extruders.

A rather recent, but promising approach is to use smoothed particle hydrodynamics (Eitzlmayr and Khinast, 2015; Robinson and Cleary, 2018; Wittek et al., 2018). The advantage of this method is that it naturally allows to compute the flow field of only partially filled extruders. However, special assumptions have to be made to recover the flow effects inside small gaps. Furtheremore, the method is still computationally expensive compared to standard methods.

Fictitious domain/immersed methods and extensions of those are also a widely spread approach. The flow solution is computed on a background mesh and the boundaries and screw surfaces are described implicitly. In (Valette et al., 2009), this method is used to compute the flow field inside a mixing section of an extruder. Standard fictitious domain methods require massive local refinement to capture flow inside small gaps, if not impossible, as shown in (Fard et al., 2012a). Therefore, they use an XFEM-based method to overcome this problem. It is used in (Fard et al., 2012b) and (Fard and Anderson, 2013) to compute distributed mixing for a wide range of screw configurations. A different approach to avoid loss of information inside the small gaps is the Body Conformal Enrichment method presented in (Ilinca and Hétu, 2010). This method has been used to compute the flow for multiple conveying and mixing screw element extruders (Hétu and Ilinca, 2013). One major drawback of fictitious domain methods is that portions of the mesh are inactive for individual time steps. Therefore, load balancing has to be used to maintain the performance for highly parallelized computations. In (Mierka et al., 2014), a method is presented where a mesh deformation technique is used to concentrate the background mesh close to boundaries, and by that, improve the drawbacks of fictitious domain methods in terms of resolving the flow in small gaps as well as load balancing.

A further approach are interface tracking methods, where boundary-conforming meshes are used. The boundary-conforming nature of this method leads to accurate results since the exact geometry is taken into account. Thus, they are extremely popular in case only steady flow is solved (Bravo et al., 2000). Computing solely the flow field of the polymer melt assuming a quasi-steady behaviour of the flow is a valid assumption. This allows to generate boundary-conforming meshes for each screw orientation of interest and solve the flow field on each mesh independently. This technique is used by (Malik et al., 2014) to investigate the effect of pressure-dependent slip boundary conditions. Furthermore, there exist several works that also include the steady-state energy equations to compute the flow and temperature field of non-Newtonian fluids inside the extruder (Ishikawa et al., 2001; Kalyon and Malik, 2007). In contrast to the flow field that can be assumed to be quasi-steady or periodic, a quasi-steady temperature field can only be seen as an approximation. The temperature solution should also take convective transport and its relation to the change of the flow domain into account. This implies that the current mesh of a given screw orientation uses the solution of the previous orientation. This requires the projection of the old solution onto the current mesh, which is computationally extremely expensive. Furthermore, it is not sufficient to consider only a few screw orientations, since the rotation between two orientations is dependent on the required time step size. Thus, meshing has to be performed for many orientations, which can be a very tedious task.

A different approach would be to update the boundary-conforming mesh to adapt to the new geometry using mesh update methods. A popular mesh update method is the Elastic Mesh Update Method (EMUM) (Tezduyar et al., 1992b; Stein et al., 2003) that treats the mesh as an elastic body. However, the scew movement in the small gaps introduces a lot of shear into the mesh, leading to mesh failure after few time steps. This even happens in case nodes are allowed to slide on the screw surface. Therefore, once again re-meshing and projections of solutions are necessary. Within the literature, there exist different methods that allow boundary-conforming meshes for multiple rotating domains, e.g., the Shear Slip Mesh Update Method (Behr and Tezduyar, 1999), sliding meshes (Bazilevs and Hughes, 2008) or multiple reference frames (Pauli et al., 2015). Unfortunately, they all rely on non-intermeshing rotating objects, which makes them unsuitable for twin-screw extruders.

Within this work, we present a mesh update method - Snapping Reference Mesh Update Method (SRMUM) - that allows to use boundary-conforming meshes for twin-screw extruders equipped with conveying elements. The mesh is updated solely based on algebraic operations making it a very efficient method. Furthermore, all advantages of boundary-conforming methods can be employed without any need for re-meshing and projection. The paper is structured as follows: First we specify the twin-screw extruder screw elements used within this work. Section 3 explains the basic idea of SRMUM. The governing equations used to describe the polymer melt inside the twin-screw extruder as well as the solution methods to solve them are presented in Section 4. Section 5 contains 2D and 3D examples validating the presented method as well as a complex temperature-dependent test case showing the applicability of the method to real applications.

2 Problem Statement: Screw Configurations

We will focus on standard conveying screw elements used inside co-rotating twin-screw extruders. The extruder consists of two screws located inside two circular barells. A 2D cross-section of a standard screw is given in Fig. 1. The two barrels overlap, and we will refer to the resulting intersections as cusp points. The shape of the barrel is defined by the centerline distance $C_{l}$ between the screw centers and the barrel radius. The barrel radius $R_{b}$ consists of the screw radius $R_{s}$ and the screw-barrel clearance $\delta_{b}$. Furthermore, we will only consider convex screw shapes.

For industrial application the screw contour is given by CAD data. However, a general definition of screw shapes is presented in (Booy, 1978). Geometric relations are derived that allow to generate a screw setup with constant screw-screw clearance $\delta_{s}$ throughout the whole rotation based on given $C_{l}$, $R_{s}$, $\delta_{b}$ and $\delta_{s}$. A slightly adapted version is presented in (Fard et al., 2012b). It will be used within this work to generate different screw designs. 3D conveying screw elements are generated by constant rotation of the 2D slices in positive axial direction. The length of a full revolution is called pitch length. One can distinguish between two types of conveying elements, namely forward- and backward-conveying, see Fig. 2. For forward-conveying elements, the 2D screw surface is rotated in the opposite direction to the screw rotation itself. In contrast, the rotation is the same for backward-conveying elements.

3 Mesh Update: Snapping Reference Mesh Update Method (SRMUM)

In this section, we will present the basic ideas of the Snapping Reference Mesh Update Method (SRMUM) with application to twin-screw extruders. As already indicated by the name, the general idea of the method is that a background mesh constantly snaps to the current geometry and by that generates a boundary-conforming mesh. In order to illustrate the functioning of the method, we will start using a single screw inside one barrel in 2D.

Due to the convex shape of the screw, one can map it onto a circle. Thus, a structured background mesh between an inner and outer circle can be used. Additionally, we discretize the screw surface and the barrel into $n_{s}$ equally sized elements. The number of elements for the screw and barrel discretization and the number of circumferential elements of the background mesh has to match. The nodes of the screw and barrel discretization are labeled using an ID ranging from 0 to $n_{s}$. In order to obtain a relation between the nodes of the background mesh to the screw discretization, every mesh point also receives an ID that is computed based on $ID=integer(\;\theta/360^{\circ}*\;n_{s}\;)$. $\theta$ is the angle enclosed by the line from the circle origin to the node. Thus, all nodes on a mesh line in radial direction have the same ID, see Fig. 3(a). In order to account for the screw rotation throughout the simulation, the IDs of the screw have to be shifted whenever the screw rotates more than $\theta=360^{\circ}/n_{s}$, whereas the IDs of the background mesh are constant, see Fig. 3(a) - 3(c).

In addition to the ID every node needs to have information about its relative position between the inner and outer circle. We denote this as $d_{rel}$. It is zero for points on the inner circle and one on the outer one. The distribution between the circles can be arbitrarily chosen, e.g., to account for boundary layer effects. Based on $d_{rel}$, the position of every node at any time can be determined by

[TABLE]

In order to extend this method to twin-screw extruders, we connect two background meshes in the intermeshing area between the two cusp points and eliminate duplicate nodes. The open question is how to update mesh points inside this area based on Eq. (1), due to the lack of an outer barrel point. We solve this by defining a middle line in this area that describes the interface between the two background meshes. The middle line will consist of points belonging to each ID of the intermeshing area replacing ${\bf x}_{barrel}(ID_{i})$ in Eq. (1). Its construction is complicated because the middle line has to ensure that the resulting mesh does not overlap or tangle. Furthermore, it has to be updated for every screw orientation. Its definition will be described in the following section.

3.1 Definition of the Middle Line for Twin-Screw Extruders

The basic idea of the middle line is that a point on that line for a certain $ID$ relates to the corresponding screw surface points ${\bf x}_{screw}(ID)$. A first, however naive, definition would be to take the average of the corresponding screw points $avg({\bf x}_{screw})$. Looking at the y-value, this results in very large elements close to the cusp points, see Fig. 4(a). Therefore, we also need to take the corresponding ID information into account. In order to simplify the notation, we will only consider the upper half of the intermeshing area. However, the lower part works exactly the same. We aim to relate the y-value not only to $avg(y_{screw})$ but also to the fictitious outer barrel $R_{barrel}*sin(CP_{rel}*\theta_{CP})$, where $CP_{rel}=ID/ID_{CP}$ is the relative distance between $\theta=0^{\circ}$ and $\theta_{CP}$. In a first approach this is done using a $y_{rel}$ defined as:

[TABLE]

Thus, $y_{midLine}$ is computed via:

[TABLE]

When comparing Fig. 4(a) and Fig. 3(b), the improvement is directly visible.

A further problem of using $avg(y_{screw})$ is that it might happen that this value is larger than the $y_{CP}$. The previous definition of $y_{rel}$ is not sufficient to resolve this problem. Thus, a $y_{rel}^{max}$ for $ID=ID_{CP}-1$ being the one closest to $ID_{CP}$ has to be found based on Algorithm 1. However, it might still happen that the y-values of IDs, $ID=ID_{CP}-1-i$, based on Eq. (1) are larger than the cusp point value. Algorithm 2 is used to compute the ID from which on $y_{rel}$ is allowed to be zero. Based on that we compute the decay of $y_{rel}^{max}$ by $\Delta y_{rel}$ so that $y_{rel}$ can be computed as $y_{rel}=y_{rel}^{max}-(ID_{CP}-ID-1)*\Delta y_{rel}$. In case $y_{rel}$ goes below zero, we relate $y_{rel}$ of Eq. (2) to $CP_{rel}$. A detailed description how this is done and the resulting $y_{midLine}$ is given in Algorithm 3.

Similar to the y-coordinate, we find difficulties using $avg(x_{screw})$ close to the cusp points, see Fig. 4(c). Depending whether the screws penetrate the left or right barrel, the middle line should follow the fictitious left or right barrel wall inside the intermeshing domain. However, it could happen that that line intersects the screws. Therefore, once again, we relate the averaged screw value and the barrel using $CP_{rel}$. A detailed description is given in Algorithm 4. Based on that a valid middle line can be computed for any orientation, see Fig. 4(d).

It is important to note, that the SRMUM is not only tailored to one specific screw configuration. It works for any screw that can be generated using Booys formula presented in Section 2.

3.2 Extension to 3D

The extension of SRMUM to 3D is rather easy. Conveying screw elements are built through rotations of 2D slices. These 2D slices can be interpreted as 2D meshes generated using SRMUM, see Fig. 5. All these 2D meshes are generated based on the same background discretization. Therefore, we can simply connect the resulting slices to build a full 3D mesh. Large rotations between two consecutive slices might lead to twisted elements. However, the element length in axial direction is restricted by the fact that it has to be small enough to cover the underlying physics correctly. For most applications this is already enough to avoid twisted elements. Additionally, it is important to note that the discretizaton of the background mesh is defined by the number of elements in screw ($n_{s}$), radial ($n_{r}$) and axial ($n_{a}$) direction, and either tetrahedral or hexahedral elements can be used.

Remark**.**

The beauty of SRMUM is that no extra equation has to be solved to update the mesh. Every node of the mesh is only updated through algebraic operations. This makes SRMUM a very efficient mesh update method. The computational time for the mesh update is a tiny fraction of the time used for the solver. Furthermore, no tedious re-meshing and expensive projection is needed.

4 Governing Equations and Solution Method

4.1 Governing Equations for the Flow and Temperature of Plastic Melt

The behavior of the molten polymer in an extruder is modeled as the flow of a viscous incompressible temperature-dependent fluid. The time-dependent computational domain, denoted by $\Omega_{t}\;\subset\;\mathbb{R}^{n_{sd}}$, is enclosed by its boundary $\Gamma_{t}$, where $t\in(0,T)$ is an instant of time and $n_{sd}$ the number of space dimensions. The velocity ${\bf u}({\bf x},t)$, pressure $p({\bf x},t)$ and temperature $T({\bf x},t)$ are governed by the incompressible Navier-Stokes equations:

[TABLE]

where $\rho$ is the fluid density and $\eta$ the dynamic viscosity. Considering a Newtonian or Generalized Newtonian fluid the following relation of the stress tensor $\bm{\sigma}$ is used to close the set of equations:

[TABLE]

The Dirichlet and Neumann boundary conditions for temperature and flow are defined as:

[TABLE]

$\left(\Gamma_{t}\right)^{i}_{g}$ and $\left(\Gamma_{t}\right)^{i}_{h}$ are complementary portions of $\Gamma_{t}^{i}$, with $i=f\mbox{ (Fluid)},\;T\mbox{ (Temperature)}$.

In order to model the shear-thinning and shear-thickening behavior of the polymer melt Generalized Newtonian models have to be used. All these models have in common, that the viscosity depends on the invariants of the rate of strain tensor $\bm{\varepsilon}$, such as the shear rate $\dot{\gamma}$:

[TABLE]

Within this work we use two different models, namely the Carreau and the Cross-WLF model.

The Carreau model is one of the most popular shear-thinning models in the plastic community (Carreau and De Kee, 1979; Bird et al., 1987). Its advantage stems from the fact that it is valid over a broad range of shear rates, particularly $\dot{\gamma}\to 0$. It is defined as:

[TABLE]

where $\eta_{0}$ is the viscosity at zero shear shear rate, $\eta_{\infty}$ is the viscosity at infinite shear rate, $\lambda$ is the relaxation time and $n$ is the power index.

The Cross-WLF model considers not only effects of shear rate but also of the temperature on the viscosity (Rudolph and Osswald, 2014). We assume that the infinite viscosity is negligible such that the Cross model can be written as:

[TABLE]

where $\tau^{*}$ is the critical shear stress at the transition from the Newtonian plateau. It is important to note, that $\eta\left(\dot{\gamma},T\right)$ depends on the temperature now. This relation is modeled via the WLF equation:

[TABLE]

with $D_{1}$ being the viscosity at a reference temperature $T_{ref}$ and $A_{1}$ and $A_{2}$ are parameters that describe the temperature dependency.

In order to relate the momentum diffusivity to thermal diffusivity of a generalized Newtonian fluid we make use of the material specific Prandtl number $Pr=\eta_{\infty}/(\kappa/(\rho c_{p}))$ (Sato et al., 2004). This is motivated by the fact that a shear-rate dependence of the thermal conductivity for non-Newtonian fluids could be shown in (Lee and Irvine Jr, 1997; Lee, 1998; Kostic and Tong, 1999). Thus, we adapt $\kappa$ based on the shear-rate dependent viscosity by keeping the Prandtl number constant throughout the simulation. We are aware that we might introduce a little too much thermal diffusivity by that due to fact that other studies showed that the Prandtl number itself is viscosity dependent (Sato et al., 2006).

4.2 Space-time Finite Element Discretization

In order to account for the transient fluid flow we need to discretize the equations in space and time. Additionally, we deal with a constantly moving and deforming domain. One natural approach that takes all that into account is the DSD/SST (Deforming Spatial Domain / Stabilized Space-Time) method (Tezduyar et al., 1992a). This formulation does not only construct the weak form of the underlying equations for the spatial- but the corresponding space-time domain. By that, we avoid the necessity of modifying the equations to account for the deforming domain.

Next, we will define the finite element function spaces for the DSD/SST method. The time interval $(0,T)$ is divided into subintervals $I_{n}=(t_{n},t_{n+1})$, where $n$ defines the time level. Setting $\Omega_{n}=\Omega_{t_{n}}$ and $\Gamma_{n}=\Gamma_{t_{n}}$, the volume enclosed by the two surfaces $\Omega_{n}$, $\Omega_{n+1}$ and the lateral surface Pn is defined as a space-time slab $Q_{n}$. Pn is described by $\Gamma_{t}$ as it traverses $I_{n}$.

For each space-time slab, the finite element function spaces for first order polynomials in space and time are defined as

[TABLE]

The stabilized space-time formulation for the incompressible Navier-stokes equations (4) and (5) without the heat equations (6) then reads as:

Given $({\bf u}^{h})^{-}_{n}$, find ${\bf u}^{h}\in(\mathcal{S})_{\bf u}^{h}$ and $p^{h}\in(\mathcal{V})_{p}^{h}$ such that:

[TABLE]

holds for all ${\bf w}^{h}\in(\mathcal{V})_{\bf u}^{h}$ and $q^{h}\in(\mathcal{V})_{p}^{h}$.

The stabilized space-time formulation for the heat equation (6) reads as:

Given $(T^{h})^{-}_{n}$ find ${T}^{h}\in(\mathcal{S})_{T}^{h}$ such that:

[TABLE]

holds for all $v^{h}\in(\mathcal{V})_{T}^{h}$.

We make use of the following notation:

[TABLE]

The stabilization parameters $\tau_{\mbox{\tiny{MOM}}}$, $\tau_{{\mbox{\tiny{CONT}}}}$ and $\tau_{\mbox{\tiny{TEMP}}}$ are based on expressions given in (Pauli and Behr, 2017). Equation (22) and (23) are solved using the Newton-Raphson method. Flow field and temperature are coupled strongly using a fixed-point interation until convergence. Note, that the stress contributions in the GLS stabilization terms (fifth term in equation (22) and fourth term in equation (23) ) are zero, since they involve second derivativ es and only first order polynomials are used. In order to improve the consistency of our method we employ a least-squares recovery technique for these terms (Jansen et al., 1999).

Remark**.**

In case of a Newtonian model the viscosity is independent of shear rate and temperature, meaning that the heat equation (6) can be decoupled from the continuity (4) and momentum (5) equation. Thus, no fixed-point iteration is necessary and a one-way coupling between flow field and temperature is used.

5 Numerical Examples

5.1 Validation Case: 2D Flow through Twin-Screw Extruder

In order to show the functioning of the newly presented SRMUM in 2D, we use a test case presented in (Fard et al., 2012a). We simulate the isothermal flow of a plastic melt through a 2D cross section of a twin-screw extruder. The screw geometry is based on Booy’s description presented in Section 2. The geometric paramters are given in Table 3. The screws rotate with $\omega_{s}=$ 60 rpm in mathematically positive direction.

The plastic melt is modeled using the Carreau model and its parameters are given in Table 3. Not considering temperature effects and taking into account that the time-scales of the momentum diffusion are very small compared to the process itself, the flow is quasi-steady or instantaneous. Therefore, we do not compare results for the whole rotation but only consider 2 screw orientations ($\theta=0^{\circ}$ and $\theta=112.5^{\circ}$). The time step size is $0.0125\;s$ or $4.5\;^{\circ}/s$. A no-slip boundary condition is used on the barrel and the rotational velocity is set as Dirichlet boundary condition on the screws. In order to show convergence of our solution we compute the flow field on 4 different meshes, see Table 3. In the following, we will refer to them by their number of elements in screw direction.

The flow behaviour is particularly complex inside the intermeshing domain (Fard et al., 2012a). Inside the small gap between the screws, their rotational velocities are in opposite direction resulting in a high pressure drop that drives the flow solution. The velocities tangential to the screws are extremely high causing high shear rates and high local differences in viscosity. It is crucial that our method computes the flow in this area correctly. Therefore, we compare the flow solution along a line in y-direction – for orientations $\theta=0^{\circ}$ at $x=-2.2754\;mm$ and for $\theta=112.5^{\circ}$ at $x=2.077\;mm$. The two lines are visualized in red in Fig. 6 for mesh 1.

The solutions of the pressure field computed on mesh 4 are shown in Fig. 7. The high pressure drop in the gap region is clearly visible. Fig. 8 shows the solution of the normalized y-velocity with $V_{p}=2\pi R_{s}\omega_{s}$ and pressure over the line in the intermeshing area for orientation $\theta=0^{\circ}$. The difference of the solutions between the coarsest and finest discretization are less than 0.5$\%$. The highest difference occurs, as expected, in the small gap where the high velocities occur. However, the structure of the solution is the same for the coarsest and finest mesh. The effect of the linear discretization can be observed for the velocity on mesh 1 in the peak region. Additionally, it is noteworthy that the results match those presented in (Fard et al., 2012a).

Fig. 9 shows the results over the line for orientation $\theta=112.5^{\circ}$. The flow pattern is more complex due to the non-symmetric geometry. Still the difference in the flow and pressure solution between finest and coarsest mesh over the entire line is rather small with roughly 1$\%$ magnitude. However, bigger differences occur in the peak velocity region. Here, the flow velocity differs by around 3$\%$. An overshoot for the pressure can be observed for mesh 1. The coarse discretization in combination with highly stretched elements and the singularity of the screw geometry cannot handle the high pressure drop. This is a special problem in 2D and is not observed in the 3D examples. An explanation could be that in 3D, the flow is not forced to go through this small gap because it can also escape into the extra dimension. Still, the coarsest solution qualitatively matches the finest one and we observe convergence with decreasing element size.

5.2 Validation Case: 3D Newtonian Flow through Twin-Screw Extruder

Within this section we aim to show that the meshes using SRMUM combined with methods presented in Section 4 enable us to compute valid 3D flow fields in co-rotating twin-screw extruders. In contrast to 2D, 3D computations are able to show the potential of the SRMUM. Due to the construction of a 3D mesh out of 2D slices it is ensured that SRMUM produces valid meshes for any orientation of a 2D slice. We compute the flow inside a single forward-conveying screw element. The mass flow through the extruder is larger than the natural mass flow of the screw element. This results in a pressure decrease from the inflow to the outflow, meaning that material is pushed, instead of transported, through the screw. This typically does not occur in industrial screw configuration, but is a common setup for experimental studies.

5.2.1 Grid Convergence Study

In a first step, we want to show convergence of our flow solution in 3D. The screw parameters are listed in Table 5. We consider corn syrup, which can be modeled as a Newtonian fluid with a viscosity of $\eta=4.7\;Pa\;s$ and density $\rho=1400kg/m^{3}$. A pressure boundary condition is used at the inflow with $p_{in}=0.2\;MPa$ and a natural boundary condition is set at the outflow. To show the mesh convergence, we use 5 different meshes where the fifth mesh is used as a reference solution, see Table 5. We compute 5 time steps with a time step size of $0.01\;s$.

As a measure of convergence we compare the averaged mass flow rate over the 5 time steps. As it is an integrated value, it gives a good estimation of the quality of the solution. The convergence plot is shown in Fig. 10. A clear convergence can be observed. It is important to note that even the solution on the coarsest mesh only differs by 3$\%$. Thus, even the coarsest solution has an acceptable error for industrial applications.

5.2.2 Validation against Experimental Results

We were able to show mesh convergence for 3D flow computations in twin-screw extruders. However, this is no proof that those results are realistic compared to experiments. Therefore, we aim to reproduce experimental data presented in (Bakalis and Karwe, 2002). A laser Doppler anemometry was used to obtain velocity measurements in a Plexiglas bench-top model of a co-rotating, self-wiping twin-screw extruder (ZSK-30, Krupp Werner $\&$ Pfleiderer, Ramsey, NJ). The screw parameters are the same as presented in Table 3 with a pitch length of 28 $mm$. Corn syrup at 30 ${}^{\circ}C$ was used as a model fluid. The density is 1400 $kg/m^{3}$ and the viscosity $4.7\;Pa\;s$. The screw speed is 120 rpm and the screws rotate in mathematically negative direction. The mass flow rate is 0.009117 $kg/s$. We model the mass flow by applying a constant velocity at the inlet of the extruder. The barrel and screws are again modeled using no-slip conditions. We use the discretization of mesh 4 from Table 5 to compute the flow field.

The region for comparison is the nip region in the x-y plane for one screw orientation, see Fig. 11. It is important to note that a different coordinate system is used for the experiments. Their x-velocity is equivalent to ours but the positive axial flow direction is in negative y-axis in the experiment and in positive z-axis in our model.

Fig. 12 shows the measured and computed x-velocity. The computed velocity magnitude range matches the experimental data. Furthermore, a similar isoline structure can be observed. The same results can be observed for the velocity in axial direction, see Fig. 13. The isolines structure shows a nice agreement. Again, the velocity magnitude range is similar, but the computational results produce slightly higher peak values. One can observe that the velocity is non-zero on the screw surface in the experiments. This could also explain the lower peak values, but contradicts the general belief of no-slip walls. This phenomenon has been discussed in (Malik et al., 2014). However, this topic will not be examined in further detail in this work.

Both the qualitative agreement with the experiment and the mesh convergence shown in the previous section indicate that the methods presented are valid and powerful tools to compute the flow in co-rotating twin-screw extruders.

5.3 Validation Case: 3D Shear-Thinning Flow through Twin-Screw Extruder

So far, we only considered Newtonian fluids. However, considering real polymer melts, it is also important to take the shear-thinning behavior into account. Hence, we conduct a convergence study using the shear-thinning fluid modeled by the Carreau model, see Table 3. To show the functioning of SRMUM for a wide range of screws we again use a different forward-conveying screw element, see Table 7. For simplicity, only a single screw element is modeled. We compare different operating conditions, ranging from a zero pressure difference between inflow and outflow resulting in maximal mass flow up to a change of the flow direction caused by a high pressure difference. We employ 4 different mesh discretization denoted by the screw discretizations – $n_{s}=180$, $n_{s}=304$, $n_{s}=500$ and $n_{s}=768$, as shown in Table 5.

Fig. 14 shows the mass flow rate $\dot{m}$ for 7 different pressure drop values ranging from 0 to 2 $MPa$. A clear mesh convergence can be observed. Looking at the volumetric distribution of the viscosity, a similar conclusion can be drawn, see Fig. 15. The results for mesh 768 and mesh 304 differ at most by 4 $\%$ for all operating conditions. This is again less than the 5$\%$ tolerance required in many engineering applications (Mierka et al., 2014). The difference between the coarsest and finest mesh are also less than 5$\%$ for pressure drop ranges from 0.0 to 1.0 $MPa$. In case the pressure drop increases further, the flow becomes more complicated since it starts to be pushed against the natural transport direction of the screw element. This is also indicated by higher shear rates and as a consequence lower average viscosities. In this situation, the coarsest discretization is not sufficiently accurate. However, in most relevant applications, the operating conditions are not set in such a way that they produce a mass flow rate going to zero.

5.4 Temperature-Dependent Flow in a 4 Element Twin-Screw Extruder

Within this section, we present solutions for the temperature-dependent flow inside a 4 element twin-screw extruder. The Cross-WLF model is used to account for the temperature-dependent shear-thinning fluid. This test case demonstrates the advantage of the presented method. The temperature solution depends mainly on convective flow transport so that solutions depend strongly on the previous time step. Therefore, it is a major advantage to use the same mesh throughout the whole computation, because no time-consuming projection of solutions onto new meshes is necessary.

We consider two different configurations – one with only forward-conveying screw elements (f-c) and a second one where the second screw element is changed to a backward-conveying element (b-c). The screw used is inspired by the ZSK 320 of the ZSK MEGAcompounder series. Inside twin-screw extruders, backward-conveying elements are used to generate a pressure build-up and furthermore, to ensure that the screw is fully filled. The screw parameters and the two configurations are given in Tables 10, 10 and 10. In order to avoid unnaturally high shear rates close to the inflow, we extend the extruder by one more screw element. The screw radius is gradually decreased until it reaches a circle at the inlet. With this, we circumvent high viscous dissipation source terms resulting in a nonphysical temperature increase in the inflow region. The same screw element is also added to the outflow. This ensures that no backflow occurs at the outlet, which would cause problems due to the zero heat flux condition we apply there. The two configurations are shown in Fig. 16.

We use a polymer that is similar to polypropylene. The parameters for the Cross-WLF model are presented in Table 12. The polymer has density $\rho=710\;kg/m^{3}$, specific heat $c_{p}=2400\;J/kg\;K$ and zero thermal conductivity $\kappa_{0}=0.18W/m\;s$. In order to compute the relation between thermal conductivity and viscosity we use a viscosity at infinite shear rate $\eta_{\infty}=7.5Pa\;s$. Furthermore, we prescribe a mass flow rate at the inlet of $\dot{m}=25560\;kg/h$. The screws rotate with $120\;rpm$ in mathematically negative direction. The time step size is $\Delta t=0.01s$. The barrell is heated with $T_{barrel}=473K+20K\;*\;z/1.5m$ and the screws are considered adiabatic. We set the initial condition for the temperature as $T_{0}=473K+15K\;*\;z/1.5m$.

In a first step, we only consider configuration 1. Screw and axial discretization of mesh 1 matches mesh 1 of Section 5.3, only $n_{s}$ is doubled to accurately represent the temperature boundary layer. We showed that for moderate pressure built-ups, this discretization yields sufficiently accurate results. In order to confirm this assumption also for the underlying example, a second mesh with double number of elements in each direction is also used, see Table 12. We compute flow and temperature fields for 20 time steps on both meshes. Fig. 17(a) shows the pressure in axial direction along the barrel for time $t=0.2s$. The temperature solution is visualized for a line in axial direction at $(x=28.9m,y=0.0m,z)$, see Fig. 17(b). This line is located in the center of the smallest gap between screw and barrel. The pressure and temperature solutions on both meshes show the same structure and the average error of temperature and pressure is less than 1$\%$. Thus, we can conclude that the discretization of mesh 1 is already sufficiently fine for the underlying problem. Therefore, only mesh 1 is used in the following to compute multiple revolutions of the screws. By that we can save a significant amount of computational time.

In a second step we want to analyze the effect of the the backward-conveying element onto the flow and temperature solution. We compute 750 time steps resulting in a total simulation time of $7.5s$. Fig. 16 and 18 show the the solution for the pressure at time $t=1.0s$. The pressure build-up effect of the backward-conveying screw is clearly visible. Instead of a positive pressure increase form inlet to outlet of $\sim 3.2\;MPa$ for config. 1, the backward-conveying element decreases the pressure drop so that we obtain a pressure decrease of $\sim 1.0\;Mpa$. Furthermore, it can be observed that the structure of the pressure solution is basically independent of the upstream screw configurations. The third and fourth screw element are the same for both configuration. The two pressure solutions match shortly after the end of the second screw element at roughly $z=0.7m$, since the outflow pressure is fixed to zero.

The evolution of the temperature inside the twin-screw extruder for both configurations is shown in Fig. 19. The rotating screws in combination with small gap sizes produces a lot of shear that increases the melt temperature through viscous dissipation. A particular increase of temperature inside the gap regions can be observed. The temperature solution reaches a periodic steady state after 13$\sim$14 revolutions. The backward-conveying element also has a strong effect on the temperature of the melt. The flow has to overcome the backflow transport property of these elements resulting in the pressure build-up and high shear which heats up the melt. This effect can be clearly observed for the given test case. In contrast to the pressure, the convective transport of the temperature means that the downstream solution is also affected. The maximum temperature for config. 1 is $7K$ lower than for config. 2. Thus, backward-conveying elements clearly increase the temperature. This result is in accordance with steady results presented in (Ishikawa et al., 2001).

6 Conclusion

We presented the Snapping Reference Mesh Update Method (SRMUM) as a mesh update method that enables us to use a boundary-conforming finite element method to compute the flow inside co-rotating twin-screw extruders. The method is based on a background mesh that constantly snaps to the actual screw configuration. The mesh update is done only by algebraic operations and no additional partial differential equation has to be solved. It circumvents tedious and time-consuming re-meshing, which makes SRMUM computationally really efficient. Furthermore, it was incorporated into a solver framework using the Deformable-Spatial-Domain/Stabilized Space-Time (DSD/SST) finite element formulation.

We validated the method for flow of a polymer melt modeled by the Carreau model inside a 2D cross section of a twin-screw extruder. 3D test cases for a single conveying element using corn syrup as model fluid were studied and mesh convergence as well as agreement with experimental data could be shown. The effect of the pressure drop on the mass flow rate and viscosity distribution was presented for a polymer melt in a conveying extruder. In this context, we were also able to obtain mesh convergence.

However, the beauty of SRMUM comes into play when looking at temperature-dependent flow inside twin-screw extruders, since only then the automatic mesh update is of great importance. Therefore, we presented flow and temperature results for a 4 element extruder using forward- and backward-conveying elements. A periodic quasi-steady state for the temperature field was obtained after several revolutions. We could show that backward-conveying elements have a much stronger heating effect on the melt than forward-conveying elements.

The automatic mesh update and continuous background discretization will also allow to investigate mixing based on advection-diffusion equations in the future. However, it is still necessary to adjust the method to allow computations for kneading elements, which is presently a subject of investigation.

Acknowledgements

The authors gratefully acknowledge the research funding which was provided by SABIC. Furthermore, we thank Christos Varsakelis for his support and ideas. The computations were conducted on computing clusters supplied by the Juelich Aachen Research Alliance (JARA) and the RWTH IT Center.