Boundary-Conforming Finite Element Methods for Twin-Screw Extruders using Spline-Based Parameterization Techniques

TL;DR

This paper introduces a spline-based meshing technique for boundary-conforming finite element simulations in twin-screw extruders, enabling efficient, accurate, and adaptable flow and temperature modeling for complex geometries.

Contribution

The paper presents a novel spline-based meshing method that generates boundary-conforming meshes efficiently at run-time, reducing memory and computational costs for complex screw geometries.

Findings

Effective boundary-conforming meshes generated for complex screw geometries.

Method reduces re-meshing costs and improves simulation accuracy.

Demonstrated applicability to 2D and 3D non-Newtonian fluid flows.

Abstract

This paper presents a novel spline-based meshing technique that allows for usage of boundary-conforming meshes for unsteady flow and temperature simulations in co-rotating twin-screw extruders. Spline-based descriptions of arbitrary screw geometries are generated using Elliptic Grid Generation. They are evaluated in a number of discrete points to yield a coarse classical mesh. The use of a special control mapping allows to fine-tune properties of the coarse mesh like orthogonality at the boundaries. The coarse mesh is used as a 'scaffolding' to generate a boundary-conforming mesh out of a fine background mesh at run-time. Storing only a coarse mesh makes the method cheap in terms of memory storage. Additionally, the adaptation at run-time is extremely cheap compared to computing the flow solution. Furthermore, this method circumvents the need for expensive re-meshing and projections of solutions making it efficient and accurate. It is incorporated into a space-time finite element framework. We present time-dependent test cases of non-Newtonian fluids in 2D and 3D for complex screw designs. They demonstrate the potential of the method also for arbitrarily complex industrial applications.

Tables

Table 1. Table 1: Carreau parameters.

${\eta }_0$	1290	$Pas$
${\eta }_{\infty }$	0	$Pas$
$n$	0.559	-
$\lambda$	0.112	$s$

Table 2. Table 3: Mesh discretization for 2D convergence study.

mesh	$n_s$ C-grid	$n_s$ separator	$n_s$ total	$n_r$	$\#$ elements	$n_{\mu }$	$n_{\nu }$ total
1	200	80	280	6	3360	12	280
2	300	160	460	10	9200	12	460
3	300	300	600	12	14400	12	600
4	300	600	900	18	32400	12	900

Table 3. Table 4: Mesh discretizations for 2D mixing elements for Config. 1.

mesh	$n_s$ C-grid	$n_s$ separator	$n_s$ total	$n_r$	$\#$ elements	$n_{\mu }$	$n_{\nu }$ total	$n_{slices}$ $\pi$	memory savings
1	150	60	210	4	1680	12	210	101	22 $\%$
2	300	120	420	8	6720	12	420	101	61 $\%$
3	600	240	840	16	26880	12	840	101	81 $\%$

Table 4. Table 5: Mesh discretization for 2D mixing elements for Config. 2.

$n_s$ C-grid	$n_s$ separator	$n_s$ total	$n_r$	$\#$ elements	$n_{\mu }$	$n_{\nu }$ total	$n_{slices}$ $\pi$	memory savings
400	140	540	8	8640	12	540	101	63 $\%$

Table 5. Table 6: Cross-WLF parameters.

$D1$	5.0e+13	$Pas$
${\tau }^{\ast }$	2500.0	$Pa$
$n$	0.29	-
$T_{ref}$	263.15	$K$
$A1$	28.32	-
$A2$	51.60	$K$

Table 6. Table 7: Geometry parameters of a 3D mixing screw element for temperature-dependent flow.

Screw radius $R_s$	0.156 $m$
Center line distance $C_l$	0.262 $m$
Screw-screw clearance ${\delta }_s$	0.004 $m$
Screw-barrel clearance ${\delta }_b$	0.004 $m$
Pitch length $p_l$	0.28 $m$

Equations

$$\Xi =\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{n+p+1}\}.$$

$$N_{i,s}(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+1}-{\xi }_i}{N}_{i,s-1}(\xi )+\frac{{\xi }_{i+s+1}-\xi }{{\xi }_{i+s+1}-{\xi }_{i+1}}{N}_{i+1,s-1}(\xi ),$$

$$\displaystyle N_{i,0}=\left\{\begin{array}[]{ll}1&\text{if }\xi_{i}\leq\xi\leq\xi_{i+1}\\ 0&\text{otherwise}\end{array}\right.,$$

$$\sum _{i=1}^n{N}_{i,p}(\xi )=1,$$

$$N_{i,p}(\xi )\ge 0,$$

$$\Xi =\left\{0,0,0,0,\frac{1}{7},\frac{2}{7},\frac{3}{7},\frac{4}{7},\frac{5}{7},\frac{6}{7},1,1,1,1\right\}.$$

$$w_{i,j}(\xi ,\eta )=N_{i,p}(\xi )M_{j,q}(\eta ).$$

$$\mathbf{x}=\sum _i\sum _j{\mathbf{c}}_{i,j}{w}_{i,j},$$

$$\mathbf{x}=\sum _{i\in {\mathcal{I}}_{\text{inner}}}{\mathbf{c}}_i{w}_i+\sum _{j\in {\mathcal{I}}_{\text{boundary}}}{\mathbf{c}}_j{w}_j,$$

$$\mathbf{x}(\xi ,\eta )$$

$$g_{22}{\mathbf{x}}_{\xi \xi }-2g_{12}{\mathbf{x}}_{\xi \eta }+g_{11}{\mathbf{x}}_{\eta \eta }$$

$$\left(\begin{array}{cc}{g}_{11}& g_{12}\\ g_{21}& g_{22}\end{array}\right)$$

$$\mathbf{u}(\xi ,\eta )$$

$$U(\mathbf{u},\mathbf{x})$$

$$\forall {\mathbf{\sigma }}_i\in [\tilde{\Sigma }{]}^2\times [{\Sigma }_0{]}^2:{\int }_{\hat{\Omega }}{\mathbf{\sigma }}_i\cdot \left(\begin{array}{c}{\mathbf{x}}_{\xi }-\mathbf{u}\\ U(\mathbf{u},\mathbf{x})\end{array}\right)\mathrm{d}\mathbf{\xi }=0,\text{s.t. }\mathbf{x}{|}_{\partial \hat{\Omega }}=\partial \Omega ,$$

$$\eta =\frac{1}{2}({\hat{\eta }}_i^l+{\hat{\eta }}_j^r)$$

$$\Sigma ={\sigma }_{\Xi }^{p_1}\times {\sigma }_{\mathcal{H}}^{p_2}\text{and}\bar{\Sigma }={\sigma }_{\bar{\Xi }}^{q_1}\times {\sigma }_{\bar{\mathcal{H}}}^{q_2}$$

$${\mathbf{\Xi }}^{p_1,p_2}={\Xi }^{p_1}\times {\mathcal{H}}^{p_2}\text{and}{\bar{\mathbf{\Xi }}}^{q_1,q_2}={\bar{\Xi }}^{q_1}\times {\bar{\mathcal{H}}}^{q_2},$$

$${\Xi }^{p_1}=\{\underset{p_1+1\text{ times}}{\underbrace{0,\dots ,0}},{\xi }_1,\dots ,{\xi }_q,\underset{p_1\text{ times}}{\underbrace{0.5,\dots ,0.5}},{\xi }_r,\dots ,{\xi }_s,\underset{p_1+1\text{ times}}{\underbrace{1,\dots ,1}}\},$$

$${\bar{\Xi }}^{p_1+1}=\{\underset{p_1+2\text{ times}}{\underbrace{0,\dots ,0}},{\xi }_1,\dots ,{\xi }_q,\underset{p_1+1\text{ times}}{\underbrace{0.5,\dots ,0.5}},{\xi }_r,\dots ,{\xi }_s,\underset{p_1+2\text{ times}}{\underbrace{1,\dots ,1}}\}.$$

$$\frac{1}{2}\underset{[0,1{]}^2}{\int }{\left(\frac{\partial ({\mathbf{x}}_i\circ {\mathbf{s}}_i)}{\partial \mu }\cdot \frac{\partial ({\mathbf{x}}_i\circ {\mathbf{s}}_i)}{\partial \nu }\right)}^2\mathrm{d}\mu \mathrm{d}\nu \to \underset{\{\widetilde{{\mathbf{c}}_i}{\}}_i}{\min },\text{subject to boundary conditions.}$$

$$\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),$$

$$\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),$$

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

$$\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)}},$$

Full text

††journal: Computer Methods in Applied Mechanics and Engineering

1 Introduction

Co-rotating twin-screw extruders are widely-used devices in the plastic-producing industry. They are used to distribute and disperse polymers and additives since they provide short residence times and extensive mixing. Typical screw elements are conveying elements that are used to transport the plastic melt from the feed section towards the die as well as kneading and mixing elements. The latter are especially tailored to account for dispersive and distributive mixing.

Performing experiments to obtain detailed information about velocity, pressure, and temperature distribution is very complex, time-consuming, and expensive, if not impossible for twin-screw extruders. This makes Computational Fluid Dynamics (CFD) using finite element analysis an appealing tool to obtain detailed information about the flow inside twin-screw extruders. However, numerical analysis of the flow is not trivial due to constantly moving domains and small gap sizes which makes meshing extremely difficult. To tackle this problem, a wide range of approaches have been proposed.

1D models have been developed in order to predict the pressure build-up [1] or fill length and specific mechanical energy [2]. Full 3D results have been obtained using commercial software like POLYFLOW which is based on a mesh superposition method that avoids re-meshing [3, 4]. A rather recent development are methods using smoothed particle hydrodynamics [5, 6]. These methods are extremely useful in case only partly filled extruders are considered.

Fictitious domain methods are also quite popular. Meshing is avoided since the actual geometry is embedded into a fixed background mesh. A classical example is the method presented in [7] that is used in [8] to compare 3D results with 1D estimates. However, in [9] it has been shown that classic fictitious domain methods lack accuracy inside the small gap regions. This drawback can be mitigated by, for example, XFEM. Other examples of fictitious domain methods are the Body Conformal Enrichment method presented in [10] or a method that tries to concentrate the background mesh at the screw interface using algebraic operations [11] . As a drawback, such methods require additional effort to describe the screw boundary exactly as well as to capture the flow effects in the small gaps. Furthermore, load balancing for highly parallel large scale computations is extremely difficult.

A different approach are boundary-conforming methods. The mesh represents the geometry exactly in terms of the underlying finite element discretization, which allows to strongly impose boundary conditions as well as the construction of high-quality boundary-layer meshes. The drawback is that the generation of a time-dependent boundary conforming mesh can be difficult. Thus, most of the simulations are based on a ”snapshot” technique. An individual mesh is generated for each screw orientation. This is a valid approach for flow computations since the flow inside the extruder can be considered to be quasi steady. This approach has been applied in many publications such as [12, 13]. Also steady temperature results have been obtained in [14, 15]. However, these methods are limited when simulating time-dependent quantities. A continuous time-dependent mesh is needed. This requires efficient mesh update techniques to avoid constant re-meshing and in turn expensive projections onto to the new mesh.

In [16], we present the Snapping Reference Mesh Update Method (SRMUM), an efficient mesh update method for twin-screw extruders. It allows the computation of time-dependent temperature results without any need for re-meshing. The method uses a structured background mesh that constantly adapts to the current screw configuration in a boundary-conforming manner. As a drawback, SRMUM is limited to convex screw geometries.

Instead, [17] intends to generate analysis-suitable spline-based parameterizations. The resulting meshes are suitable for numerical analysis using Isogemetric Analayis (IGA) as introduced in [18]. The approach is based on Elliptic Grid Generation, which can handle nonconvex screw geometries. The method is PDE-based, which makes it suitable for gradient-based shape optimization by simply adding the parameterization PDE-problem as a constraint to the shape optimization formulation, which can then be differentiated with respect to the design variables [19]. Furthermore, classical finite-element meshes have been extracted from the spline parameterization and used to compute first flow results inside idealized twin-screw compressors [20].

In this work, we aim to adapt the general concept of spline-based meshing developed for twin-screw compressors to twin-screw extruders and combine it with certain aspects of SRMUM. This allows us to generate flexible, high quality, time-dependent and FEM-suitable meshes.

The approach is based on taking a large number of cross sections through a volumetric extruder geometry leading to $2$D extruder contours at a large number of rotational angles $\theta_{i}$ (see Figure 1). A spline parameterization $\mathbf{x}_{i}$ is generated for the $2$D profiles at every $\theta_{i}$ using Elliptic Grid Generation techniques.

The spline parameterizations are evaluated in a large number of points to generate a mesh ‘scaffolding’ to adapt a SRMUM-like background mesh to the current screw shape via interpolation. A volumetric mesh is then built from the planar meshes by stacking them in the $z$-direction.

Once a spline parameterization has been generated, it is relatively cheap to generate finer meshes, e.g. for refinement studies. Furthermore, function evaluations in the spline mappings $\mathbf{x}_{i}$ can be tuned using a ’control mapping’ in order to account for special properties of the mesh, like orthogonality at the screw surface. Since only a ‘scaffolding’ is extracted, it is not necessary to save the whole $3$D mesh for all possible time instances which would be very time- and memory-consuming. Furthermore, the computational cost of the mesh update at each time step during the simulation is only a fraction of the time spent for solving the flow solution.

The paper is structured as follows: Section 2 explains the meshing concept. This includes the generation of a spline parameterization up to the analysis-suitable mesh. The underlying physics and equations required to compute the flow inside twin-screw extruders are stated in Section 3, as well as the space-time finite element solution method. Section 4 contains $2$D and $3$D examples of flows through twin-screw extruders that aim to validate as well as highlight the advantages of the presented method.

2 Grid Generation

Since the geometry pipeline provides no more than a description of the twin-screw profile (typically in terms of a point cloud), the first step towards numerical simulation is the generation of a planar computational mesh. For the mesh generation, we employ the following two-step approach: at every discrete angle $\theta$, we first generate a spline-based geometry description which we then evaluate it in a large number of discrete points with known connectivity in order to yield a classical mesh. The advantages of a spline-based description is the possibility to fine-tune the properties of the mesh by appropriately choosing the evaluation abscissae, which follow from a precomputed control mapping.

In the following, we give a brief introduction to spline-surfaces after which we discuss the basics of the numerical algorithm employed in the meshing process and its applications to the target geometry.

2.1 B-Splines

As described in [21], B-splines are piecewise-polynomial functions that can be constructed so as to satisfy various continuity properties at the places where the polynomial segments connect. Their properties are determined by the entries of the so-called knot vector

[TABLE]

The knot vector is a nondecreasing sequence of parametric values $\xi_{i}\subset[0,1]$ that determine the boundaries of the segments on which the spline-basis is polynomial. Selecting some polynomial order $p$, the $p$-th order spline-functions $N_{i,p}$ are constructed recursively, utilizing the relation (with $\tfrac{0}{0}\equiv 0$)

[TABLE]

starting from

[TABLE]

and iterating until $s=p$. The support of basis function $N_{i,p}$ is given by the interval $\mathcal{I}_{i,p}=[\xi_{i},\xi_{i+p+1}]$ and the amount of continuous derivatives across knot $\xi_{j}$ is given by $p-m_{j}$, where $m_{j}$ is the multiplicity of $\xi_{j}$ in $\mathcal{I}_{i,p}$. In practice, $\xi_{1}=0$ is repeated $p+1$ times as well as $\xi_{n+p+1}$ such that $\xi_{1}=\ldots=\xi_{p+1}=0$ and $\xi_{n+1}=\ldots=\xi_{n+p+1}=1$. As a result, the corresponding basis $\sigma=\{N_{1,p},\ldots,N_{n,p}\}$ forms a non-negative partition of unity on the entire parametric domain $[0,1]$, that is:

[TABLE]

with

[TABLE]

for all spline functions $N_{i,p}$ [18]. Figure 2 shows the $p=3$ B-spline basis resulting from the knot-vector

[TABLE]

The extension to bivariate spline bases is now straight-forward: given two univariate bases $\sigma_{\Xi}=\{N_{1,p},\ldots,N_{n,p}\}$ and $\sigma_{\mathcal{H}}=\{M_{1,q},\ldots M_{m,q}\}$, we build a bivariate basis $\Sigma=\{w_{i,j}\}_{(i,j)\in\{1,\ldots,n\}\times\{1,\ldots,m\}}$, by means of a tensor-product, where

[TABLE]

The values contained in $\Xi$ and $\mathcal{H}$ (without knot-repetitions) hereby become the boundaries of the polynomial segments, which can be regarded as the counterparts of classical elements.

We construct the mapping of a B-spline surface as follows:

[TABLE]

where the $\mathbf{c}_{i,j}\in\mathbb{R}^{2}$ are referred to as the control points. We refer to the $\mathbf{c}_{i,j}$ with $i\in\{1,n\}$ or $j\in\{1,m\}$ as the boundary control points while the remaining $\mathbf{c}_{i,j}$ are called inner control points. Bluntly put, the local density of spline basis functions is determined by the local number of knots. As such, a cleverly chosen knot-vector can be utilized to properly resolve important features of the geometry. Furthermore, as $\mathbf{x}$ is a linear combination of the $w_{i,j}\in\Sigma$, it will inherit the local continuity properties of the basis. Therefore, many geometrical features can be better captured by a clever choice of the knot multiplicities in $\Xi$ and $\mathcal{H}$.

2.2 Spline-Based Meshing Techniques

In the following, we shall drop the tensor-product indexing used in (10) and replace it by a single global index. As such, the mapping operator will be of the form

[TABLE]

where $\mathcal{I}_{\text{inner}}$ and $\mathcal{I}_{\text{boundary}}$ refer to the index-set of inner and boundary basis functions, respectively. The mapping $\mathbf{x}$ is a function from the parameter domain $\hat{\Omega}=[0,1]^{2}$ to the physical geometry $\Omega$, whose boundaries $\partial\Omega$ follow from a collocation of the employed spline-basis $\Sigma=\{w_{1},\ldots,w_{N}\}$, resulting from the tensor product knot vector $\bm{\Xi}$, to the input point cloud.

For the proper selection of $\bm{\Xi}$ and the subsequent collocation, we utilize the approach illustrated in [19]. Our methodology consists of a stabilized least-squares fit of the points against the basis and an adaptive reselection of $\bm{\Xi}$ based on the local magnitude of the projection residual. The process is repeated until the residual is deemed sufficiently small. Upon completion, the restriction of $\mathbf{x}$ to $\partial\hat{\Omega}$ parameterizes $\partial\Omega$. Note that the above procedure essentially selects the boundary control points in (11) such that $\mathbf{x}|_{\partial\hat{\Omega}}$ is a proper collocation of the input point cloud (hence the word boundary control points). However, the inner control points are unknown at this stage. As such, the task of any parameterization algorithm is to properly select the $\mathbf{c}_{i}$ in (11), such that the resulting parameterization is folding-free.

A computationally inexpensive, yet often sufficiently powerful method is transfinite interpolation [22]. The inner control points follow from an interpolation of the southern and northern boundaries in combination with an interpolation from east to west. Let $\bm{\gamma}_{e},\bm{\gamma}_{w},\bm{\gamma}_{s}$ and $\bm{\gamma}_{n}$ parameterize the four segments of $\partial\Omega$. In a spline-based setting, the mapping operator is constructed as follows:

[TABLE]

where the $\mathbf{p}_{i,j}$ denote the corners at $\mathbf{x}(i,j)$. The symbolic parameterization from (12) constitutes a recipe for determining the $\mathbf{c}_{i}$ in (11). They can be computationally inexpensively determined by a $L_{2}$-projection or a collocation over the Greville-abscissae [23] corresponding to $\bm{\Xi}$. Whenever $\mathbf{x}$ is an O-type parameterization, we simply remove the nonexistent pair of boundary contours ($\bm{\gamma}_{n},\bm{\gamma}_{s}$ or $\bm{\gamma}_{w},\bm{\gamma}_{e}$) from (12), as well as the $\mathbf{p}_{i,j}$ and perform a unidirectional interpolation.

Transfinite interpolation is an $\mathcal{O}(N)$ operation but does not guarantee a folding-free mapping. In order to handle the complex characteristics of twin-screw extruders, it is desirable to have a more powerful parameterization technique in our arsenal.

A second class of parameterization methods are PDE-based, notably approaches based on the principles of Elliptic Grid Generation (EGG). The main purpose of EGG is to compute a mapping $\mathbf{x}:\hat{\Omega}\rightarrow\Omega$ such that the components of $\mathbf{x}^{-1}:\Omega\rightarrow\hat{\Omega}$ are harmonic. This is accomplished by imposing the following system of PDEs on $\mathbf{x}$:

[TABLE]

where the $g_{ij}$ denote the entries of the metric tensor corresponding to $\mathbf{x}$. We have

[TABLE]

A justification of this approach is based on the observation that if $\partial\Omega$ satisfies certain regularity conditions, $\mathbf{x}$ will be bijective (and hence, folding-free) [24]. As such, any sufficiently accurate approximation of $\mathbf{x}$ will be bijective as well.

2.3 Numerical Implementation

A basic numerical algorithm to approximately solve (13) with globally $C^{\geq 1}$-continuous spline basis functions is discussed in [17]. Here, we follow a similar approach with auxiliary variables, first introduced in [25], that allows for $C^{0}$-continuities in $\xi$-direction (whose purpose shall become apparent in subsection 2.4). We approximately solve (13) with a Galerkin approach, introducing the auxiliary variable $\mathbf{u}$ that serves the purpose of reducing the continuity requirements of the basis in $\xi$-direction from $C^{1}$ to $C^{0}$.

Let

[TABLE]

where the $\tilde{w}_{i}\in\tilde{\Sigma}$ are basis functions taken from the auxiliary spline basis $\tilde{\Sigma}$.

In order to discretize (13), we express the $\xi$-derivatives in terms of $\mathbf{u}$ and introduce a scaling. Consider

[TABLE]

with $\epsilon>0$ a small term that ensures numerical stability. In order to determine the $\mathbf{c}_{i}$ in (11) and $\mathbf{d}_{i}$ in (15), we solve the system

[TABLE]

where $\Sigma_{0}=\{w_{i}\in\Sigma\enskip|\enskip i\in\mathcal{I}_{\text{inner}}\}$. For a memory-efficient algorithm to tackle the root-finding problem (17), we refer to the aforementioned publication [25].

Upon solving (17), we are in the possession of $\mathbf{x}$ and $\mathbf{u}$. Here, $\mathbf{x}$ parameterizes $\Omega$ and $\mathbf{u}$ serves no further purpose and can be discarded. Unlike the exact solution of (13), its numerical approximation may fold due to numerical inaccuracies. Since a defect is a direct cause of insufficient numerical accuracy, it can be repaired by increasing the local density of spline functions wherever the defect is located and recomputing a mapping operator from the enriched spline basis. In practice, we have observed folding to be a rare occurrence.

2.4 Applications to the Geometry

For the parameterization of the main geometry, we choose the topology illustrated in Figure 3. The topology consists of two $C$-type parameterizations connected by a third parameterization (blue), henceforth referred to as the ’separator’. The topology is designed in such a way that upon rotation, the rotor lobes slide along the grid in the parametric domain while the grid is held fixed at the casings. Furthermore, the casings are held arc-length parameterized while the parametric properties of the rotor lobes depend on the rotation angle.

As such, the first step in building a three-patch parameterization is finding a valid $O$-grid parameterization of a single rotor plus casing. We take the circular casing and map it onto a straight line segment. A valid parameterization is one in which the vertical isolines do not intersect under the mapping (see Figure 4).

To find a valid lobe parameterization, we employ the matching algorithm from [19] to the rotor and casing point clouds. The points are matched based on their Euclidean distance in a hierarchical fashion and matched points receive the same parametric value. Upon completion, we have a monotone reparameterization function $\xi^{\prime}:[0,1]\rightarrow[0,1]$ and we assign the parametric value $\xi_{i}=\xi^{\prime}(\hat{\xi}_{i})$ to the $i$-th point in the rotor point cloud. Here, the $\hat{\xi}_{i}$ result from a chord-length parameterization of the point cloud. The steps described above are carried out for $\theta=0$ only. The reparameterization function can be reused for different rotational angles by performing an appropriate periodic shift that is based on the value of $\theta$.

The coordinates of the points along with their parametric values are utilized for a collocation of the left and right rotor $O$-grids to a suitable spline basis. Here, the knot vector is chosen such that the collocation residual is below a user-defined threshold (see [19]). Upon completion, the boundary curves are known and we utilize interpolation (see (12)) to form folding-free $O$-grid parameterizations for both rotor casing pairs. Next, we cut the left and right rotor $O$-grids at the parametric values that correspond to the casing CUSP points to form the left and right rotor $C$-grids.

We repeat above steps for a large number of discrete angles $\Theta=\{\theta_{1},\ldots,\theta_{n}\}$ taken from the interval that corresponds to a rotation after which the initial position is again obtained (typically the interval $[0,\pi]$). Upon completion, we have a database that is filled with $C$-grid parameterizations for all discrete angles $\Theta$. We combine the $C$-grid boundaries along with the rotor lobe parts that were cut from the rotor-casing $O$-grids in order to generate one boundary description of the separator, as shown in Figure 3. Clearly, a spline collocation of the northern and southern boundaries has to be performed with a locally $C^{0}$-continuous spline basis, due to the spiked nature of the input point clouds. Therefore, we assign the parametric value $\xi=0.5$ to the CUSP-point and utilize a knot vector with $p$-fold internal knot repetition at $\xi=0.5$ for the $\xi$-direction.

The next step is to assign suitable parametric values to the points contained in the eastern and western input point clouds. Unfortunately, we have observed a chord-length parameterization to lead to unsatisfactory results. In order to ensure that similar parametric values are assumed on either sides of small gaps in the geometry, we again apply the matching algorithm from [19]. In contrast to the $C$-grid case, we do not keep one side chord-length parameterized but let both sides float. If two points $p_{i}^{l}$ and $p_{j}^{r}$ have been matched, we assign the parametric value

[TABLE]

to both of them, where the values of $\hat{\eta}_{i}^{l}$ and $\hat{\eta}_{j}^{r}$ correspond to chord-length parameterizations. Upon completion, we have two reparameterization functions that we utilize to assign parametric values to the input point clouds, similar to the $C$-grid case.

We perform the above steps for each $k$-th discrete angle $\theta_{i}$. Typically we take $k=5$. A global reparameterization function is constructed by blending each $k$-th reparameterization function over the $\theta$-interval. Blending serves the purpose of achieving a degree of smoothness in the parametric properties of the separator as a function of $\theta$, which would be lost if we reparameterized at every $\theta_{i}$, due to the discrete nature of the reparameterization algorithm.

Upon completion, we start filling our database with separator parameterizations utilizing the EGG algorithm from (17). The auxiliary spline space $\bar{\Sigma}$ is of Raviart-Thomas type [26] on each of the macro elements [27] $\bm{\xi}\in[0,0.5]\times[0,1]$ and $\bm{\xi}\in[0.5,1]\times[0,1]$ with $C^{0}$ interface coupling. To be precise, let

[TABLE]

be built from the knot vectors

[TABLE]

respectively. Here, the $p_{i}$ and $q_{i}$ denote the polynomial orders used in each direction. If

[TABLE]

we take $\bm{\bar{\Xi}}^{q_{1},q_{2}}=\bar{\Xi}^{p_{1}+1}\times\mathcal{H}^{p_{2}}$, with

[TABLE]

In our examples, $\Sigma$ is bicubic, i.e., $p_{1}=p_{2}=3$.

The root-finding problem is initialized with a transfinite mapping (see (12)). Convergence can be further accelerated using multigrid techniques (see [17]). We fill our database in a hierarchical fashion. As soon as the database contains separator parameterizations for a sufficiently large subset of angles from $\Theta$, we utilize interpolation in order to generate initial guesses for the remaining angles to further speed up the process.

After the algorithm has finished, we have a three-patch parameterization for the geometry at every $\theta_{i}$. A classical mesh can be generated by evaluating the spline mapping in a large number of uniformly-spaced points. The properties of the mesh can be further tuned by using a control mapping. We have observed this to be especially important for the mesh quality in the separator. By the tuple $(\mathbf{x}_{i},\mathbf{s}_{i})$, we denote the mapping operator and the corresponding control mapping for the separator at $\theta=\theta_{i}$. The control mapping $\mathbf{s}_{i}(\mu,\nu)=\sum_{i}\tilde{\mathbf{c}}_{i}\tilde{w}_{i}(\mu,\nu)$ is a function from the unit square onto itself. We choose its control points by solving the following minimization problem:

[TABLE]

The cost function from (22) maximizes orthogonality in the composite mapping $\mathbf{x}_{i}\circ\mathbf{s}_{i}$. The integral is approximated with a discrete Riemann sum. The boundary conditions follow from the requirement that $\mathbf{s}_{i}$ be a mapping from the unit square onto itself. We only let $\mathbf{s}_{i}$ slide in $\eta$-direction (see Figure 5) as this leads to better results for tube-like shaped geometries such as the separator. A sufficient linear condition for the aforementioned requirement is easily formulated and imposed as a constraint on (22). The minimization problem is tackled with an SLSQP [28] algorithm while the gradient of the cost function is approximated with finite differences. To speed up convergence, control mappings of previous slices can be utilized to initialize the next control mapping. We compute a control mapping for the separator at all discrete angles $\theta_{i}\in\Theta$. We do not compute control mappings for the $C$-grids (i.e., the control mapping is the identity).

Upon completion, we build a mesh ’scaffolding’ for the separator from the composite mapping $\mathbf{x}_{i}\circ\mathbf{s}_{i}$ by evaluating it in a large number of uniformly-spaced points. By $n_{\mu,\nu}$, we denote the number of elements in $\mu$ and $\nu$ direction, respectively. An adaptive background mesh is generated by adding additional vertices via linear interpolation within the mesh scaffolding point cloud. The vertex positions are computed on the fly and do therefore not have to be stored explicitly (see Figure 6). For the background mesh, we denote the number of elements in radial and screw direction by $n_{r}$ and $n_{s}$ with $n_{r}\geq n_{\mu}$ and $n_{s}\geq n_{\nu}$, respectively. The ratios $n_{r}/n_{\mu}\geq 1$ and $n_{s}/n_{\nu}\geq 1$ are hence based on a trade off between memory requirements and grid quality, where ratios closer to $1$ typically lead to better meshes (for a given number of elements) while increasing the memory requirements. The converse holds for larger ratios, whereby it is important to take $n_{\mu,\nu}$ sufficiently large so that the resulting mesh ’scaffolding’ is sufficiently close to the folding-free $\mathbf{x}_{i}\circ\mathbf{s}_{i}$ and therefore does not fold itself.

For the C-grids, the mesh scaffolding cannot fold due to an insufficient number of elements in $\mu$-direction, thanks to the precomputed reparameterization function (see Figure 4). Hence, we take $n_{\mu}=1$ and add vertices until $n_{r}(C)=n_{r}(\text{separator})/2$ to yield a boundary-conforming mesh on both sides of the CUSP points. The storage requirements for the C-grids are thus fully determined by $n_{\nu}$, whose value is chosen to properly resolve the rotor boundaries. Typically, for the C-grids, we simply take $n_{s}=n_{\nu}$.

In order to construct a 3D mesh, we linearly connect the 2D meshes, where $n_{a}$ determines the number of elements in $z$-direction. Every slice in $z$-direction knows its initial angle, which is used to interpolate between the stacked planar slices of the ‘scaffolding’.

3 Governing Equations and Solution Method

3.1 Governing Equations for Flow and Temperature of Plastic Melt

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

[TABLE]

where $\rho$ is the fluid density. The stress tensor $\bm{\sigma}$ is used to close the set of equations:

[TABLE]

with $\eta$ being the dynamic viscosity. It is constant for a Newtonian fluid and varies for Generalized Newtonian models with respect to temperature $T$ and shear rate $\dot{\gamma}$.

Using single-phase Navier-Stokes equations implies the assumption that the extruder is fully filled. However, in practice it might occur that the flow channels inside the twin-screw extruder are only partially filled. In order to account for this, one would have to extend the model to also include the air phase, which results in solving a multi-phase flow problem. Popular numerical methods for multi-phase flow are the levelset [29] or the volume of fluid method [30]. For simplicity, only single-phase flow will be considered within this work.

The temperature $T({\bf x},t)$ inside the extruder is governed by the heat equation:

[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)}$.

It is of course true that plastic melts exhibit viscoelastic behavior. We assume, however, that this is negligible in this context, since the residence time in the twin-screw extruder is small. We use Generalized Newtonian models to account for shear-thinning and shear-thickening behavior of the polymer melt. Using Generalized Newtonian models implies 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 a very popular shear-thinning model in the plastic community [31]. It is defined as:

[TABLE]

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

The Cross-WLF model considers along with the effects of the shear rate also the influence of temperature on the viscosity [32]. We neglect the infinite viscosity such that the Cross model is defined as:

[TABLE]

where $\tau^{*}$ is the critical shear stress at the transition from the Newtonian plateau. $\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.

Remark**.**

In most of the available literature, Stokes equations are used to model the fluid inside a twin-screw extruder. However, Reynolds numbers around 0.1 occur in our applications s.t. the Stokes equations are not necessarily a valid assumption any more. Thus, we model the fluid using the Navier-Stokes equations which also include the advective term. It is noteworthy that the nonlinearity introduced due to the advective term is small compared to the one introduced by using Generalized Newtonian models.

3.2 Space-time Finite Element Discretization

We need to discretize the equations in space and time in order to capture the transient fluid flow. Additionally, the rotating screws introduce a constantly moving and deforming domain. A natural approach accounting for all of these requirements is the DSD/SST (Deforming Spatial Domain / Stabilized Space-Time) method [33]. In DSD/SST, the weak form is written not only over the spatial domanain, but instead the space-time domain. Thus, we do not need to modify the equations to account for the deforming domain.

In the following, 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. We set $\Omega_{n}=\Omega_{t_{n}}$ and $\Gamma_{n}=\Gamma_{t_{n}}$. Thus, a space-time slab $Q_{n}$ is defined as the volume enclosed by the two surfaces $\Omega_{n}$, $\Omega_{n+1}$ and the lateral surface Pn. Pn is described by $\Gamma_{t}$ as it traverses $I_{n}$.

We use first-order interpolation for all degrees of freedom. Thus, a SUPG/PSPG stabilization technique is used in order to fulfill the LBB condition [34]. The finite-element interpolation and weighting function spaces for velocity, pressure and temperature for every space-time slab are defined as

[TABLE]

The stabilized space-time formulation for equations (23) and (24) then reads:

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

[TABLE]

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

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

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

[TABLE]

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

The following notation is used:

[TABLE]

The stabilization parameters $\tau_{\mbox{\tiny{MOM}}}$, $\tau_{{\mbox{\tiny{CONT}}}}$ and $\tau_{\mbox{\tiny{TEMP}}}$ are based on expressions given in [35]. We use a Newton-Raphson method to solve equation (41) and (42). The flow and temperature fields are coupled strongly using a fixed-point iteration until convergence is reached. The stress contributions in the SUPG/PSPG stabilization terms (fifth term in equation (41) and fourth term in equation (42) ) are zero, since they involve second-order derivatives. We improve the consistency of our method by employing a least-squares recovery technique for these terms [36]. In order to solve the resulting linear system of equations within each Newton iteration, we use a GMRES solver with an ILUT preconditioner.

4 Numerical Examples

4.1 2D Isothermal Flow

We want to show that simulations based on the new method for grid generation produces the same high-quality results as already established methods like SRMUM or XFEM. Therefore, we simulate the isothermal flow of a plastic melt, described by the Carreau model, in a 2D cross section of a twin-screw extruder. The Carreau parameters are given in Table 2. The screw geometry is generated based on an adapted version of Booy’s description as presented in [9, 16]. The screw parameters are given in Table 2. The rotation speed of the screws is $\omega_{s}=60$ rpm in mathematically positive direction. We set the rotational velocity on the screws as a Dirichlet boundary condition and a no-slip condition on the barrel. In [9], Stokes flow was used. Since we use Navier-Stokes equations to model the flow, we employ a density of $\rho=1\;kg/m^{3}$ in order to make the results comparable.

We employ a time step size of $\Delta t=0.00625s$ or $2.25\;^{\circ}/s$. The resulting flow can be considered as quasi-steady or instantaneous because the time-scales of the momentum diffusion are very small compared to the process itself. The resulting Reynolds number inside the small gap region is $Re=0.000003$. Therefore, it is sufficient to compare the solution at individual screw orientations. Two screw orientations, namely $\theta=0^{\circ}$ and $\theta=112.5^{\circ}$ have been used in [16]. It is demonstrated that $\theta=112.5^{\circ}$ is the more complex orientation. Thus, we will only compare results for that orientation.

In order to verify the general applicability of the method, we use 4 different mesh resolutions. One advantage of the presented method is that we have to construct the spline parameterization only once and can then extract the ’scaffolding’ for all desired mesh resolutions. The region of special interest is the intermeshing area between the two screws. A high pressure drop drives the flow solution resulting in high velocities and high shear rates, see Figure 7(b). Thus, resolving this part correctly is extremely important. Therefore, we only refine our mesh inside this region and keep the mesh resolution inside the two C-grids constant. The only exception is the first mesh, where we slightly reduce the resolution for the C-grid in order to avoid highly stretched elements. We evaluate the spline-parameterization of the C-grids as well as the separator for all mesh points on the screw, meaning $n_{s}=n_{\mu}$. Inside the separator we use a ’scaffoling’ with 12 elements in $\mu$-direction. This number is kept constant throughout the refinement of the background mesh. All mesh quantities are given in Table 3. In the following, we will refer to them by the total number of elements in screw direction ’$n_{s}$ total’.

We compare the velocity results inside the intermeshing area along a line in $y$-direction at $x=2.077\;mm$. The line is visualized in red for mesh 280 in Figure 7(a). As a reference solution, we use flow results computed with the same time step size on a SRMUM mesh with 1000 elements in screw and 20 elements in radial direction. The plot over line for the normalized y-velocity $y_{vel}/V_{p}$, with $V_{p}=2\pi R_{s}\omega_{s}$ is shown in Figure 8. Looking at the overall plot over line in Figure 8, we can not observe major differences between the results computed on different meshes compared to the reference solution. The velocities differ by less than $1\;\%$. However, the peak velocities differ more. Figure 8 shows that this difference decreases with increasing mesh resolution. The results are also in accordance with the results presented in [9].

4.2 2D Temperature-Dependent Flow inside Mixing Elements

In the previous section, we have shown that the meshes based on the spline-based parameterization technique produce valid results. However, the true advantage over SRMUM only becomes apparent in the context of non-convex screw shapes - shapes that SRMUM cannot handle at all. We use two different screw configurations that are inspired by the screw design that has already been given in Table 2. The resulting screw geometries are given in Figure 9.

Config. 2 can be considered slightly more complex since the screw geometry has sharper edges compared to Config. 1. We employ several meshes for Config. 1 that will be later used in order to perform a mesh convergence study. Mesh details are given in Table 4 and 5, including the memory savings obtained by only storing a ’scaffolding’. Note, that the mesh for Config. 2 is slightly finer than the second mesh of Config. 1. This is due to the more complex shape of the screw profile, which necessitates using more grid points to properly capture the sharpe edges.

Similar to the previous section, we evaluate the spline parameterization along with the computed control mapping. Also in this non-convex case, the reparameterization is used to ensure that the C-grid parameterization with transfinite interpolation does not fold. Inside the separator, we again evaluate the spline control mapping composition in a uniform grid comprised of 12 elements in $\mu$-direction and $n_{\nu}=n_{s}$ in $\nu$-direction. The resulting interpolated fine meshes inside the separator are given for different angles for Config. 1 in Figure 10 and Config. 2 in Figure 11. In order to evaluate the mesh quality, we use the scaled Jacobian determinant calculated by the ’Mesh quality’ filter in paraview [37]. We obtain highly distorted elements in particular in the vicinity of concave corners, but the scaled Jacobian determinant is never negative. This is extremely important because the finite-element discretization can cope with highly distorted elements but fails in case of negative determinants.

In the following, we aim to demonstrate the advantage of the presented meshing method. In case of isothermal flow, it is not crucial to have matching discretization for consecutive screw orientations due to the quasi-steady behavior of the flow. However, taking temperature effects into account results in equations where a solution at the next time step strongly depends on the previous one. Therefore, we simulate a temperature-dependent plastic melt inside the two 2D mixing elements. The melt is modeled using the Cross-WLF model. The parameters are given in Table 6. Note, that the parameters have been selected for testing purposes but are in inspired by those of polypropylene. The plastic melt has density $\rho=700\;kg/m^{3}$, specific heat $c_{p}=2400\;J/(kg\;K)$ and thermal conductivity $\kappa=10.0\;W/(m\;s)$. The screws rotate in mathematically positive direction with $\omega=60$ rpm. For the flow, we set a no-slip condition on the barrel and the rotational velocity as Dirichlet condition on the screws.

Concerning the temperature, the screws are considered to be adiabatic. On the barrel we set a Dirichlet temperature condition as $T_{barrel}=473K+(10x)/0.03K/m$ and by that, generate a linear increase of the temperature from the left to the right barrel. The initial condition is $T_{0}=473K$. It is important to note that this constitutes a test case that has no connection to a real twin-screw extruder application. It merely demonstrates the validity of the method, as cold/hot melt is constantly pushed from left to right and vice versa. We use different time step sizes to show that the solutions are independent of the time step. Thus, we have to interpolate linearly between slices for most of the resulting screw orientation. Assuming that each individual slice is bijective, a linear interpolation in $z$-direction will lead to a bijective volumetric mesh.

Config. 1:

We investigate the temperature field for Config. 1 for different time steps on different meshes. In a first step we neglect the effects of viscous dissipation to solely demonstrate how the hot and cold melt is pushed from left to right and vice versa. We use three different meshes to show that the results are mesh independent. Therefore, we compare the temperature results along a line between the lower and upper cusp point at time $t=1.625s$. The time step size is $\Delta t=0.00625s$. The results are given in Fig. 12(a). Mesh convergence is clearly visible. The results on the coarsest mesh already show very good agreement to the ones obtained on the finest mesh. In the following, we investigate the choice of an appropriate time step size. We compute the temperature results on mesh 2 for four different time steps, 12(b). We can observe a difference between results for the different time step sizes especially inside the small gap region. However, even for $\Delta t=0.0125s$ we obtain results with a reasonable error for industrial applications. The difference between results obtained with $\Delta t=0.0625s$ and the finest time step size $\Delta t=0.0015625s$ are less than $1\;\%$.

Figure 13 shows the temperature field at different angles computed on mesh 2 with $\Delta t=0.003125s$. As already described, cold melt is transported by the screws from left to right. Inside the recessed portion of the screw in the lower part of the left barrel, cold melt is transported into the intermeshing area, see Figure 13(a). At time $t=1.675s$, the cold fluid has been transported into the intermeshing area and starts to be convected into the lower part of the right barrel. The flow direction inside the small gaps between the screw is negative, pushing warmer melt into this area. In the following time steps, $t=1.725s$ and $t=1.775s$, more and more warm melt is pushed through the small gap into the lower part of the intermeshing area, leading to a temperature increase there. The temperature field emerges smoothly in time showing the good quality of the underlying meshes as well as numerical methods.

In the following, we include viscous dissipation. Thus, we expect the melt to heat up due to high shear rates inside the small gaps. Again, we analyze the effect of the mesh and time step size. The results are shown in Fig. 14. Similar to the case without viscous dissipation, we see a good convergence in time and space. However, the difference of the solution between the coarsest and finest time step increases especially inside the small gap region. This shows the additional complexity introduced to the simulation by viscous dissipation. The temperature field at different angles computed on mesh 2 with time step $\Delta t=0.003125s$ is shown in Figure 15. Comparing the results with and without viscous dissipation, we can clearly observe the effects of the heating of the melt due to viscous dissipation. At time $t=1.675s$, heated melt inside the small gap is pushed into the lower part of the intermeshing area. This continues for the following time steps, leading to a major increase of temperature in that area. The temperatures are more than four degrees above those resulting from neglecting viscous dissipation. However, the convective behavior of the melt is similar to the case with no viscous dissipation.

Config. 2:

For Config. 2, we analyze the temperature and flow behavior at an earlier stage, see Figure 16. Based on the numerical experiments for Config. 1, we use a time step of $\Delta t=0.003125s$. In this configuration, the cold melt is not transported inside the recessed portion, but is pushed into the intermeshing area due to a larger gap region between screw and barrel, see Figure 16(a). Between time $t=1.55\;s$ and $t=1.6s$, the melt is pushed upwards out of the cold temperature part, which decreases the temperature in the upper part between the screws. At time $t=1.675s$ we can observe an increase of temperature in the lower part of the intermeshing region. The direction of flow inside the screw-screw gap has changed sign and the melt is pushed downwards through the gap. Again, the melt is heated up especially inside this small gap due to high shear rates. All these observations are in line with what one would expect by purely looking at flow results.

The results for the two configurations show the potential of the presented spline-based meshing approach. We are able to generate high-quality meshes for extremely challenging moving domains. Only generating a certain amount of slices and interpolating all other meshes seems to be a valid approach. Due to the space-time approach of the method, this unsteady 2D test case is already a proof of concept for 3D. Furthermore, it is noteworthy that the time spent for updating the mesh was less than $0.1\;\%$. The simulations have been run on 24 cores on the Intel Xeon based RWTH cluster using an MPI parallelization. For example, computing one revolution for Config. 1 on mesh 2 using a time step of $\Delta t=0.00625s$ takes $115s$. The time spent for the mesh update was only $0.02s$.

Remark**.**

The appeal of the presented methodology becomes apparent when performing mesh refinement studies as well as time step size studies. The spline parameterization for the screw only has to be generated once at the beginning. For a mesh refinement study only the ’scaffolding’ needs to be regenerated which is a simple evaluation of the spline parameterization. A time step refinement study is even simpler. Given a ’scaffolding’ for instances of $\Theta$ in the interval $[0,\pi]$, we can simply compute the mesh at any time instance by interpolating between the generated instances. Thus, it is not necessary to generate any new mesh in case one aims to adapt the time step size.

4.3 3D Application Case

Within this section, we aim to show the functioning of the spline-based parameterization technique for real 3D applications. Once again, we consider the temperature-dependent flow of a plastic melt through a complex screw geometry. The plastic melt is modeled using the Cross-WLF model. The model parameters employed are based on a polypropylene from the product portfolio of a leading raw material manufacturer. The parameters are given in Table 9. The plastic melt has density $\rho=710\;kg/m^{3}$, specific heat $c_{p}=2400\;J/(kg\;K)$ and thermal conductivity $\kappa_{0}=0.5\;W/(m\;s)$.

The 2D screw geometry cross section is a simplified combination of the two 2D configurations of the previous section, see Figure 17(a). The screw geometry parameters are given in Table 9. The 3D setup is shown in Figure 17(b). For simplicity we only consider a single screw element. The screws rotate in mathematically positive direction with $\omega=120$ rpm. The computational domain is extended at the inflow and outflow. By that, a solely positive velocity in $z$-direction is guaranteed at the outflow. This is important since we set a natural boundary condition and thus, negative velocities might cause instabilities. Furthermore, this setting is also similar to a industry-like twin-screw extruder where the screw section leads to the die. The extended inflow circumvents high shear rates in the inflow area, which would result in an unnaturally high temperature increase. The extension is achieved by relaxing the 2D screw surface over a distance of $0.28\;m$ to a circle with a radius of $0.06\;m$. For the mesh, we simply interpolate between the original 2D mesh and a structured mesh between the two circles. We use a background mesh with $n_{s}=220$ elements on the screw, split into 70 elements inside the separator and 150 for the C-grid. $n_{r}=12$ elements are used in radial direction. A full pitch length is discretized using 200 elements in axial direction and the inflow and outflow extensions are discretized with 50 elements each. The ’scaffolding’ generated by the spline parameterization consists of all points in circumferential direction, $n_{s}=n_{\nu}$ and $n_{\mu}=16$ inside the separator. In order to account for the rotation, the ’scaffolding’ is evaluated for 101 equally distributed instances of $\Theta$ in the interval $[0,\pi]$. Storing only the ’scaffolding’ instead of the full mesh for each instance results in memory savings for the mesh of $75\;\%$. All mesh parameters are listed in Table 9.

We simulate a scenario with a mass flow rate of $\dot{m}=25560\;kg/h$. This is achieved by setting a uniform inflow velocity at the inlet. The barrel is heated with $T_{barrel}=473K+10K\;*\;z/1.12m$ and the screws are considered to be adiabatic. The inflow temperature is $T_{inflow}=473\;K$. In order to obtain a good initial condition for the temperature, we compute a steady solution where we neglect the viscous dissipation term and increase the conductivity by a factor of 10. The time step size is $\Delta t=0.005s$ which is chosen based on the results of the previous section. The results have been computed on 360 cores on the Intel Xeon based Juelich cluster. Computing one revolution takes $5832s$, whereas the part spent for the mesh update is only $0.77s$ which shows the efficiency of the presented meshing approach.

The flow results are presented in Figure 18. The screw is able to transport a higher mass flow than the one set at the inflow. This results in an increase of pressure over the screw element, see Figure 18(b). The velocity in $z$-direction is given in a plane at $y=0.074\;m$ in Figure 18(a). A strong backflow behavior can be observed in the small gaps between the screws. This is desirable for twin-screw extruders since it increases their mixing capabilities. Inside the extruder we observe Reynolds numbers up to $Re=0.1$, which justifies using Navier-Stokes equations. The temperature field reaches a periodic state after roughly 10 revolutions. The temperature results for half a rotation are shown in Figure 19. The melt is heated inside the high shear regions in the small gaps between the screws. Furthermore, we can observe how cold melt is pushed forward inside the high velocity regions, especially inside the wide parts of the intermeshing region.

5 Conclusion

Within this work, we presented an efficient and robust meshing strategy that allows to use boundary-conforming finite element methods to compute the unsteady flow of plastic melt inside co-rotating twin-screw extruders. It is suitable for arbitrarily-shaped screw geometries. The method is a combination of spline-based meshing techniques based on EGG and SRMUM. 2D spline-based geometry descriptions are generated for a certain number of screw orientations. The spline description is evaluated at discrete points to obtain a point cloud that is used as a ’scaffolding’ to adapt a structured background mesh to the actual screw configuration. In contrast to algebraic grid generation, the advantage of the proposed approach is two-fold: on the one hand, it allows for fine tuning of the mesh properties by employing a control mapping (see Section 2.4). On the other hand, refinement studies are easily accomplished by evaluating the composite spline-mapping in a increasing number of points. Additionally, storing only a limited number of points instead of a full 3D mesh rotation saves a lot of memory. The actual mesh update at run-time is very cheap and requires only a fraction of the time spent for the actual solve of the flow solution.

A finite element method based on the Deformable-Spatial-Domain/Stabilized Space-Time (DSD/SST) finite element formulation was used. A 2D test case simulating the flow of an isothermal polymer in a twin-screw extruder cross section served as a validation case for the presented framework including the new meshing strategy. Convergence of the solution has been demonstrated as well as accordance to results from literature. Furthermore, two complex, mixing-element-like, screw shapes have been used to show the robustness of the meshing technique in 2D. We computed unsteady temperature-dependent flow of the plastic melt inside the aforementioned screw shapes. The use of a space-time finite element method already proofs, that the presented method is capable of computing unsteady flow results using boundary-conforming meshes in 3D. We additionally computed the flow of a temperature-dependent plastic melt in 3D for a complex screw shape in a single element twin-screw extruder to further demonstrate the great potential of the presented approach.

In the future, we aim to compute flow results for a larger variety of screw designs. This will be used in order to compare the quantities like residence time distributions or mixing behavior. The advantages of the method already shown for unsteady temperature results can also be exploited by computing solutions of advection-diffusion equations, in order to characterize mixing behavior of the individual screw design.

Acknowledgements

The authors gratefully acknowledge the research funding which was partly provided by the MOTOR project that has received funding from the European Unions Horizon 2020 research and innovation programm under grant agreement No 678727. The computations were conducted on computing clusters supplied by the Juelich Aachen Research Alliance (JARA) and the RWTH IT Center.