Combining Boundary-Conforming Finite Element Meshes on Moving Domains Using a Sliding Mesh Approach

For most finite element simulations, boundary-conforming meshes have significant advantages in terms of accuracy or efficiency. This is particularly true for complex domains. However, with increased complexity of the domain, generating a boundary-conforming mesh becomes more difficult and time consuming. One might therefore decide to resort to an approach where individual boundary-conforming meshes are pieced together in a modular fashion to form a larger domain. This paper presents a stabilized finite element formulation for fluid and temperature equations on sliding meshes. It couples the solution fields of multiple subdomains whose boundaries slide along each other on common interfaces. Thus, the method allows to use highly tuned boundary-conforming meshes for each subdomain that are only coupled at the overlapping boundary interfaces. In contrast to standard overlapping or fictitious domain methods the coupling is broken down to few interfaces with reduced geometric dimension. The formulation consists of the following key ingredients: the coupling of the solution fields on the overlapping surfaces is imposed weakly using a stabilized version of Nitsche's method. It ensures mass and energy conservation at the common interfaces. Additionally, we allow to impose weak Dirichlet boundary conditions at the non-overlapping parts of the interfaces. We present a detailed numerical study for the resulting stabilized formulation. It shows optimal convergence behavior for both Newtonian and generalized Newtonian material models. Simulations of flow of plastic melt inside single-screw as well as twin-screw extruders demonstrate the applicability of the method to complex and relevant industrial applications.


Introduction
In this paper we present a stabilized finite element formulation for fluid and temperature equations that allows to couple boundary-conforming discretizations of individual moving domains at common interfaces. In the world of finite element analysis for flow problems, the representation of the computational domain plays an important role, especially for complex moving domains in 3D. The numerical solution requires high-quality meshes to ensure good approximation properties on the one hand, and a proper geometric resolution of the given domain on the other hand. The question of how to balance these two aims has led to a variety of approaches and methods. Broadly speaking, these can be placed into one of two general categories: (1) boundary-conforming meshes that are aligned with the domain boundary and (2) unfitted, fictitious or overlapping methods, which all describe techniques in which the actual domain is embedded into a static background mesh or individual meshes arbitrarily overlap. In the former case, a mesh has to be generated based on given geometrical data. This can be an expensive task both in terms of human resource -since manual intervention is often necessary -as well as in terms of computer resources. At the same time, this approach ensures full control over the resolution and can easily include expert knowledge, e.g., it is possible to generate high-quality boundary layers or to refine special regions of interest already in advance. Furthermore, the imposition of Dirichlet boundary conditions is straightforward. In case of moving domains, e.g., for fluid-structure interaction the mesh can be updated in order to adapt to the moving boundary. What may seem like a drawback is, at the current state of the art, easily covered using mesh update methods based on radial basis functions [1,2] or the Elastic Mesh Update Method (EMUM) [3,4]. A more advanced extension of EMUM based on fiber-reinforced hyperelasticity is pesented in [5]. In addition, there exist a broad range of specialized mesh update methods for specific applications such as the Shear Slip Mesh Update Method (SSMUM) for rotating components [6,7] or its extension to large translation, the Virtual Ring Shear Slip Update Method (VRSSMUM) [8]. For rotating screw machines, methods that automatically adapt to the moving screw domain have been developed in [9,10,11]. Yet, even the most sophisticated method has its limitations: They occur when boundary deformations are too large or result in topological changes. At this point, the only standard option that is available is continuous remeshing during the simulation; the price to pay is the extensive use of computational resources along with a loss of accuracy of the simulation due to the necessary mapping between the individual meshes. With regard to unfitted methods, a variety of options such as level-set methods [12] or immersed boundary methods [13] are available. In the following, we will focus on methods using overlapping meshes, as these are relevant to the further development of the paper. For overlapping mesh methods, meshes are generated around static or moving objects inside the domain -e.g., the rotating object -and are then embedded into a static background mesh that covers the whole computational domain. They have been introduced under the name of Chimera [14] or Overset methods [15]. The solutions on the domains are coupled weakly using a Dirichlet/Neumann coupling [16]. The key advantage of simple mesh generation is somewhat offset by the disadvantage of the resulting in additional coupling steps. A general overview of these methods can be found in [17]. Fictitious domain methods using Lagrange multipliers that enrich the finite element function space (XFEM) are also wide spread, and have been applied to numerous complex applications [? 18, 19, 20, 21]. Instead of enriching the functions space, Nitsche's method [22] has been used in [23] to enforce the interface coupling as well as boundary conditions weakly. These Nitsche-based methods form the basis of cutFEM, since two meshes intersect or are cut by each other, see, e.g., in [24,25]. Based on this, a multimesh approach has been developed that allows to couple arbitrary many overlapping meshes [26,27]. Noteworthy are also formulations of the finite cell method that in addition use Nitsche's method [28,29]. In terms of rotating screw machines, standard fictitious domain methods accompanied with a mesh deformation technique have been developed in [30]. All the fictitious domain methods have in common that they require the computation of the cut between either the background mesh with the underlying geometry or the overlapping meshes. This is not trivial in 3D. Furthermore, load balancing is challenging for highly parallel large scale computations. Bazilevs and Hughes present an approach that combines the best from both worlds for cases including rotating with prescribed rotation [31]. An individual, boundary-conforming mesh is generated both around the rotating object and for the outer domain. In contrast to cutFEM, both domains do not overlap but share a common interface. The two meshes slide over this interface throughout the rotation, resulting in a non-matching interface one dimension lower than the actual domain. The two solution fields are coupled using Nitsche's method similar to cutFEM. The method has been successfully applied to large-scale simulations of wind turbines [32]. Note, however, that this method does not handle interfaces that change over time.
The method that will be presented in this work follows Bazilevs and Hughes [31] in the sense that it allows for boundary-conforming meshes wherever possible and couples those meshes at common interfaces. As before, these interfaces may change throughout the simulation when the meshes slide along the interface. Like Bazilevs and Hughes, we only consider overlapping meshes at boundary interfaces -e.g., a 2D cross section in 3D -thus strongly simplifying the computation of cut cells and avoiding degenerated polyhedrons. What is new is that we also take into account non-matching interfaces caused by rigid structures that exist only one side of the common interface. Here, we make use of weak imposition of boundary conditions in order to match the idea of the coupling condition [33,34]. Furthermore, we have expanded the formulation to generalized Newtonian material models and to conjugate heat transfer. The method is especially suited for applications, where individual complex moving domains are modularly stacked together in various orders. In this paper, we consider the example of extruders. The paper is structured as follows: In Section 2, we give a more detailed motivation for the presented method as well as the problem setting of the sliding mesh approach. Moreover, the governing equations as well as the coupling conditions are formulated. Section 3 states the weak stabilized formulation for each domain. Additionally, a detailed description of the weak imposition of the coupling conditions and boundary conditions using Nitsche's method is presented. This is followed by a discussion of the practical implementation. Numerical examples are presented in Section 4. Spatial convergence analysis for viscous and convective flow regimes for the 2D Taylor-Green flow, as well as for a complex 3D flow of plastic melt through a twin-screw extruder kneading element with two discs, are conducted. They show the optimal convergence behavior of the presented method. In Section 5, we apply our method to compute the temperature-dependent flow inside rotating screw machines. These test cases demonstrate the high potential of the method for complex moving time-dependent simulations of realistic industrial applications. Finally, we draw conclusions in Section 6.

A Sliding Mesh Approach for Non-Isothermal Fluid Flow
The main motivation for the present work are situations where the domain of the simulation consists of individual parts that share a common interface. Furthermore, the individual parts can rotate or move such that the domain boundaries slide over each other at the common interface. Very often, the individual domains show complex geometric features. In case of flow simulation, it is extremely important to accurately represent these geometric features, e.g., small gaps, with a proper discretization. Examples are single or twin-screw extruders (SSE/TSE) that consist of individual screw parts stacked behind each other. The resulting individual domains are characterized by extremely small gaps between screw and outer barrel. In order to capture the flow effects -especially large gradients in pressure -correctly, it is of utter importance to resolve the small gaps. Thus, boundary-conforming meshes can be beneficial. Special boundary-conforming mesh update techniques have been developed for twin-screw extruders in [9,10].  Fig. 1 shows the example of a so-called kneading element as it can be found in a twin-screw extruder. This kneading element consists of two fixed discs that are staggered by 90 • , see Fig. 1(a). Structured boundary-conforming meshes can be designed for the individual discs, shown in green and red. However, they don't match at the common interface Γ 1,2 f f . Fig. 1(b) shows the two individual meshes at the common interface as well as the overlapping part in blue. Within the next sections, we will in detail describe a method that enables us to couple the individual domains at non-matching interfaces.

Domain Coupling Using a Sliding Interface
We consider a setup in which the overall time-dependent fluid domain Ω t is divided into individual spatial subdomains Ω i t , i = 1, ..., n d , where n d is the number of subdomains, see Fig. 2. It holds that Ω t = n d i=1 Ω i t and i, j f f are denoted by Γ i SI and Γ j SI . Each subdomain Ω i t is discretized using a finite element mesh T i . We use a classical continous finite element approximation space for each subdomain K i is a tetrahedral or hexahedral element of T i . Q p denotes the polynomials of order p in each direction on K i . In principle the polynomial order can vary on each subdomain. However, only linear polynomials Q 1 will be used within this work.

Governing Equations and Coupling Conditions
We model the fluid as a viscous, incompressible fluid on the moving domain Ω t ⊂ R n sd , with n sd being the spatial dimension. As already mentioned, Ω t consists of n d disjoint subdomains Ω i t which are enclosed by their boundaries Γ i t , where t ∈ (0, T) is an instance of time. The velocity u i , pressure p i and temperature T i in every subdomain Ω i t are governed by the incompressible Navier-Stokes and heat equations in convective form: where ρ is the fluid density, b the gravity vector, κ the thermal conductivity and c p the specific heat capacity. ∂(·) ∂t represents the time derivative in the arbitrary Lagrangian Eulerian (ALE) frame and u i ALE is the domain mesh velocity. A detailed derivation of the ALE description can be found in [35,36]. The Cauchy stress tensor σ is defined as: with ε(u i ) being the strain-rate tensor and η the dynamic viscosity. η is constant for Newtonian fluids, and for Generalized Newtonian models varies with respect to temperature T and shear rateγ.
The latter is defined asγ Within this work we use two different models, namely the Carreau and the Cross model with WLF correction.
The Carreau model [37] states: where λ is the relaxation time, n is the power index, η 0 is the viscosity at zero shear rate and η ∞ is the viscosity at infinite shear rate.
The Cross model [38] is: where τ * is the critical shear stress at the transition from the Newtonian plateau. Furthermore, we also want to model the influence of temperature on the viscosity. This is done by making the viscosity at zero shear rate η 0 dependent on temperature via the WLF correction: where D 1 is the viscosity at a reference temperature T re f and A 1 and A 2 are parameters that describe the temperature dependency.
The Dirichlet and Neumann boundary conditions for flow and temperature fields are: , where superscripts f and t denote flow and temperature, respectively, and Furthermore, we need to conserve mass, momentum and energy over the common interface Γ i, j t f f between two subdomains Ω i t and Ω j t . The mass conservation in combination with a noslip condition requires continuity of velocity and temperature: Momentum and energy conservation are obtain by demanding equal surface tractions and heat fluxes: · denotes the jump in a function over the interface (Γ t ) i, j f f and is defined as · = (·) i − (·) j .

Finite Element Discretization
In the following, we discretize the equations (2) -(4). We use a SUPG/PSPG type stabilized formulation for each subdomain Ω i t following [39]. It can also be interpreted as a variational multiscale formulation as presented in [40]. The weak form is derived by multiplying the continuity equation (2), momentum equation (3) and heat equation (4) with test functions q i , w i and v i respectively. Next, we choose appropriate finite dimensional test and trial functions (1) and integrate by parts taking the boundary conditions into account. We define the L 2 -inner product on Ω i The stabilized weak formulation for the temperature equation is where I f w h , q h ; u h , p h and I t v h ; T h are the interface consistency or jump flux terms. Note that we did not include the enforcement of the interface coupling condition yet. This will be done in the next section. The first line of equations (18) and (21) are the standard Galerkin terms. The remaining lines of equations (18) and (21) are the residual-based stabilization terms, where the residuals r h CON , r h MOM and r h TEMP are defined as: The stress contributions in equations (24) and (25) involve second order derivatives. In case only first order polynomials are used, we recover these terms by using a least-squares recovery technique [41]. This improves the consistency of the stabilized method especially for highly viscous problems. The terms in line three of equation (18) as well as of equation (21) are the SUPG stabilization terms. They are used to stabilize the formulations for convection-dominated problems. The first term in line four of equation (18) adds artificial diffusion to stabilize the continuity equation. The second term is the PSPG term which is needed to make the formulation inf-sub stable, since we use equal order polynomial for all degrees of freedom. The stabilization parameters τ MOM , τ CONT and τ TEMP are based on expressions given in [39,42].

Nitsche Coupling at Common Interface
Following [25], we can now formulate the Nitsche coupling. First, we have to define the weighted interface average operator: k i and k j are real positive weights. Typical choices are one-sided weightings, meaning k i = 1, k j = 0, or an averaged weighting also denoted as {·} m where k i = k j = 0.5 as used, e.g., in [31]. However, this might lead to an unbalanced weighting between strains in case generalized Newtonian models are used. This is due to possible large jumps in viscosity between elements. Thus, we will use a weighting that balances the difference in viscosity, given as: . We use the same weighting for the heat equation but substitute η by κ. In case of constant value for η or κ, it results in the standard averaged weighting.
Using the relations given in (26), we can reformulate the interface jump terms (19) and (22) under the assumption that and by incorporating the flux condition (16) to: Using a Nitsche coupling inspired by [25,31] we can extend the interface jump term for the fluid I f w h , q h ; u h , p h (see Eq. (28)), and obtain the following stabilized formulation for the flow coupling condition at the common interface Γ f f : The terms in the second line are the so-called adjoint-consistency terms. They fulfill the coupling condition (15) and ensure mass conservation over the interface. Note that the terms involving the pressure trial (first line) and test function (second line) are skew-symmetric. This ensures that the form is stabilityneutral without violating the adjoint-consistency of the formulation. The second to last line is a consistent penalization term that ensures coercivity of the formulation and furthermore, ensures mass conservation over the interface. We define τ f SI as: with ∂ξ ∂x is the inverse Jacobian of the element mapping between reference and physical domain and α is a stabilization parameter. The last line in equation (30) is an upwinding stabilization that controls instabilities when mass is transported from one subdomain into another one [25].
The coupling for the temperature at Γ f f follows the same Nitsche approach already applied to the flow. Thus, similarly to Eq. (30), we extend I t v h ; T h ( see Eq. (29)) and obtain: We define τ t SI similar to (31) as:

Weak Imposition of Boundary Conditions Using a Multimesh Technique
Within this section, we will discuss the case in which the interface between two domains does not fully match geometrically, i.e., Γ SI Γ i f f ∩ Γ j f f . In this case, we only want the fluid to flow through the overlap, see Fig. 2. In the remaining part of the boundary interface we apply a no-slip condition to impose the movement of the underlying rigid parts of the other subdomain. Thus, we need to be able to apply a Dirichlet boundary condition in the part of Γ SI that is not connected to the other subdomain, namely Γ SI \Γ i, j f f . The Dirichlet boundary conditions will be imposed weakly also using Nitsche's method [22,33]. This is done by adding the following terms to the stabilized weak form of the flow equations: The same concept is also applied to weakly impose Dirichlet temperature boundary conditions on Γ SI \Γ i, j f f :

Algorithmic Details
In order to evaluate integrals on Γ i,j f f as they appear in (30) and (33) it is necessary to compute the cuts between boundary elements K i on Γ i SI and K j on Γ j SI in order to define quadrature rules. In contrast to embedded mesh techniques where the dimension of the cut cells is equal to the dimension of the domain, the computation of cut cells for our approach is one dimension lower than that of the computational domain. Thus, we only need to compute cut cells in 2D. The computation of the cut cells is performed using CGAL [43]. First, we detect all possible collisions of boundary elements of the subdomains at all common sliding interfaces using bounding boxes. The candidate cut elements are then checked for collision using geometric predicates. In the following, we compute the intersection of two elements K i and K j . K i ∩ K j is convex such that we can simply triangulate the cut polygon to define standard Gaussian quadrature rules on the resulting sub-triangles. The open question is how to determine the quadrature rule for the weak imposition of Dirichlet boundary conditions on Γ SI \Γ i, j f f . Defining quadrature rules for all elements on the boundary that are not cut is trivial. However, we also have to take the partly cut elements into account. The uncut part of an element is not necessarily convex, which makes sub-triangulation non-trivial (see Fig. 3). Thus, we will make use of the multimesh concept of nested quadrature introduced in [26]. Considering Fig. 3, we can interpret K i as an element on Γ i SI and K j as an element on Γ j SI such that K i ∩ K j is the resulting cut element. We aim to integrate on the uncut part of K i denoted as K i \K j . As described before, we are able to define quadrature rules on K i ∩ K j since we can present it as a set of triangles. The idea of Johansson et al [26] is to not explicitly define a quadrature rule for K i \K j . Instead they use the inclusion-exclusion principle of combinatorics: to derive Thus, we can simply integrate on K i \K j by integrating on K i and also on K i ∩ K j but using negative weights. Extending this concept to the integration of weak Dirichlet conditions on the sliding interface, we can simply integrate the terms on Γ SI \Γ i, j f f by integrating on Γ SI using standard quadrature rules and then simply integrate the terms on all cut elements on Γ f f using negative weights.
Remark. Summing up and neglecting terms everywhere on the boundary due to the negative quadrature weights sounds like a computationally expensive process at first. However, the terms of the weak imposition of boundary condition (17) and those for the Nitsche coupling at the common interface (30) have a very similar structure. Thus, basically no extra computational cost due to the subtracting arises during assembly, since the terms are already computed for the Nitsche coupling.

2D Taylor-Green Flow
We study the convergence behavior of our method computing the 2D Taylor-Green problem [44]   This problem has already been used to test approaches based on Nitsche's method for embedded as well as cut mesh methods [25,27,34]. We divide the domain into four equally sized square domains (see Fig. 4(a)).
The two subdomains on the upper left and lower right are discretized with a mesh with n elements in each direction, and the off-diagonal ones with a mesh with m elements in each direction. For the convergence study, we start with n = 4 and m = 3 and refine seven times. In order to compare our results, we also use one mesh that is not subdivided. The initial number of elements for the reference mesh in one direction is 8. We set the analytical solution as Dirichlet boundary condition on the outer boundary. In order to have a constant pressure level, we set the analytical solution for the pressure in the lower left corner. For the viscous case with η = 0.1, we compute a steady solution. It is obtained by setting the body force components b to the negative time derivative of the analytical flow solution. The steady solution for pressure and velocity is shown in Fig. 4. For the convection-dominated case, we compute 10 time steps with time step size ∆t = 0.00025. We discretize the time derivative using a Backwards Differentiation Formula of first order (BDF1) [45]. The very small time step is used in order to prevent that time discretization errors affect the convergence study. As initial condition we set the analytical flow solution (Eq. (39)) for t = 0. We analyze the L 2 error norm on the whole domain as well as the interface coupling error norm The error norms for the velocity should converge with second order and the pressure errors norms with first order [46,24].

Spatial Convergence Study
In a first step, we consider the convergence behavior for convection-dominated flow. The convergence rates are shown in Fig. 5. For the L 2 error on the whole domain we compare the results with those on a matching, uncut mesh. Furthermore, we use two different Nitsche stabilization parameters α = 10 and α = 30. These values are inspired by numerical studies performed in [25,34]. We observe that the domain as well as interface error norms converge with the expected or even higher convergence rates. The higher rates can be explained by the high regularity of the solution. Furthermore, Fig. 5(c) shows that we obtain optimal convergence for the mass conservation across the interface, which indicates that the upwinding scheme as well as the scaling of our stabilization parameters works as expected.
We observe a large difference in the error for the matching and the non-matching case, see Fig. 5(a) and 5(b). This can be explained by the fact that we plot the norm over the minimal element length that is the same for the matching and non-matching case. However, the initial number of elements for the non-matching case in two subdomains is only 3 instead of 4 (refer to Fig. 4(a)). Thus, the overall number of degrees of freedom is always smaller, which explains the offest in the error norm. We also analyze the spatial convergence for the viscous case. Fig. 6 shows the convergence behavior in all error norms. Similar to the convective case, we observe optimal convergence rates for the domain as well as interface error norms. Having a closer look at the error in mass conservation over the interface (see Fig.  6(c)), we observe that we obtain slightly smaller errors for α = 30 compared to α = 10. Thus, α = 30 will be used for the following test cases. These optimal spatial convergence rates indicate the correct scaling of our stabilization parameters for different flow regimes. We have not tested the coupling for the temperature. However, the structure of the heat equation is the same as for the momentum equation, so that similar convergence results can be expected.

3D Flow in the Kneading Element of a Twin-Screw Extruder
In this section, we extend the spatial convergence study to a more complex flow problem in 3D. Specifically, we analyze the flow inside a twin-screw extruder kneading element with two discs of length 20 mm that are staggered by 90 • . The design of the underlying screw is based on Booy's law [47,48] and the parameters are given in Table 1.        The two individual kneading elements are separately discretized using the Snapping Reference Mesh Update Method (SRMUM) as presented in [9]. Thus, we obtain two non-matching surface meshes at the interface between the two individual kneading elements that will be coupled using the presented sliding mesh approach. The computational domain is the same as used for the motivation in Section 2, see Fig.  1(a). The resulting interface is shown in Fig. 1(b). We use the Carreau model in order to account for the non-Newtonian behavior of the plastic melt inside the extruder. The Carreau parameters are given in Table 2. We compute the steady flow solution for a screw rotational speed of ω s = 60 rpm and a pressure difference between inlet and outlet of ∆p = 0.5MPa. A no-slip condition is used on the barrel, and we set the rotational speed as Dirichlet condition on the screw surface. For the spatial convergence study, we use five different meshes. The mesh parameters are given in Table 3: n s , n r , and n a describe the number of elements on the screw, in radial, and in axial direction, respectively [9]. For this test case, we also have to include the weak imposition of boundary conditions at the interface. The movement of the underlying screw from the other kneading element has to be accurately described on both sides of the interface Γ SI \Γ f f . The solution for the velocity component in z-direction on two planes computed on mesh 3 is shown in Fig. 7. We observe that the velocity contours are smooth across the non-matching interface. Furthermore, we can compare the velocity magnitude on the two sliding interfaces Γ i SI , see Fig. 8. We would like to show that the velocity of the opposite screw movement is set correctly. Looking at the velocity magnitude contours, we observe circular contour lines in those areas that are aligned with the opposite screw body. This is in perfect agreement with the fact that screw velocity magnitude increases linearly in radial direction of the screw center. Thus, we can conclude that the weak imposition of Dirichlet boundary condition using the multimesh technique works correctly. Finally, we want to have a closer look at the interface error norm, which we compute based on Eq. (41). The results for the spatial convergence analysis conducted on five consecutively refined meshes (see Table 3)    are shown in Fig. 9. Similar to the results shown in Section 4.1.1, we obtain optimal convergence rates for pressure and velocity. For this sliding mesh setup in 3D using a generalized Newtonian fluid model, this validates the weighting proposed in Section 3.2 as well as the scaling for the Nitsche stabilization parameters. This test case shows very good results, so that we can use the proposed method for more complex application cases.

Application Cases
In the following, we will apply the sliding mesh approach to two relevant physical applications in the plastics manufacturing industry. The first one is the computation of the temperature-dependent flow of a plastic melt inside a twin-screw extruder with several kneading blocks. The second one considers the flow inside single-screw extruders with varying design.

Temperature-Dependent Flow of Plastic Melt in Twin-Screw Extruder
We consider the temperature-dependent flow of plastic melt through a twin-screw extruder section with different screw elements. We will simulate the temperature evolution inside the extruder over several revolutions starting from a constant initial condition. The non-Newtonian behavior of the plastic melt is modeled by the Cross model with WLF correction, which allows to take the temperature effects into account. The model parameters are chosen based on a polypropylene from a portfolio of a raw material manufacturer, see Table 5. The melt has density ρ = 710 kg/m 3 , specific heat c p = 2900 J/kg K, and thermal conductivity κ 0 = 0.2W/m s. Similar to [9], we use the Prandtl number Pr = η ∞ c p /κ 0 to relate the momentum diffusivity to thermal diffusivity [49]; the Prandtl number used is Pr = 145000.   The twin-screw extruder is built up by seven individual screw sections -artificial relaxation sections at the beginning and the end with a length of 28 mm, where the screw shape is decreased quadratically to a circular shape, two forward-conveying elements with pitch length p l = 28 mm, one conveying element with p l = 20 mm, one kneading element 45/5/25 (45 • staggering angle, 5 discs and 25 mm length) and one kneading element 60/2/16. The screw centers are at x = ±13.1 mm and y = 0.0 mm, respectively and the inflow plane is at z = 0 mm. The 2D screw shape is again based on Booys' law, see Table 4, and the overall screw setup is shown in Fig. 10. All individual screw elements including each disc of a kneading element are discretized using SRMUM. This allows to update the mesh at each time step in an efficient way without any need for re-meshing. The number of elements in screw direction and radial direction are n s = 216 and n r = 10. The overall number of elements in axial direction is 634. This has proven to be an appropriate discretization in mesh studies presented in [9,10] for similar simulations. The individual mesh building blocks are coupled at the common interfaces using the presented sliding mesh approach. We simulate a mass flow rate ofṁ = 23.88kg/h which is achieved by setting a uniform Dirichlet inflow condition. A no-slip condition is set at the barrel and a natural boundary condition at the outflow. The screws rotate in mathematically positive direction with ω s = 120 rpm. The resulting streamlines at t = 6s are shown in Fig. 10. For the temperature, we set a uniform inflow temperate at T in = 473K. The screws are considered adiabatic and the barrel is heated using a Dirichlet boundary condition with T barrel = 473K + 5/z max · 5K. Furthermore, we set the initial temperature of the plastic melt to T 0 = 480K. The time step size is set to ∆t = 0.005s based on numerical studies presented in [9,10]. The temperature distribution inside the extruder reaches a quasi-steady periodic state after 10 revolutions.  (e) t = 5.62s (f) t = 6.0s Figure 12: Temperature distribution in the xz-plane at y = 6mm inside the twin-screw extruder with different screw sections. Fig. 11 shows the temperature distribution between the two discs of the kneading element 60/2/16 at t = 6.54s. We can observe higher temperatures especially inside the small gap regions due to viscous dissipation. Furthermore, at the common interface the temperature distributions match well. This indicates the validity of the presented sliding mesh approach also for the heat equation. Fig. 12 shows the temperature distribution over time in an xz-plane at y = 6mm. At t = 0.62s and t = 1.0s the influence of the initial condition is still visible. However, cold melt is already pushed into the extruder and a reduction of the temperature due to barrel heating and cooling can be observed. Additionally, a temperature increase in the vicinity of the small gaps as well as kneading elements due to viscous heating is visible. The overall temperature decreases over time compared to the initial condition. The temperature distribution after t = 1s is already close to the quasi-steady state. The biggest differences occur in the outflow region. It is noteworthy, that the average temperature in any cross section is higher than the given temperature at the corresponding barrel position. This clearly shows the importance of viscous heating effects in the extruder and is also illustrated in Fig. 11. The computed temperature distributions are also qualitatively in accordance to steady-state results for different kneading elements presented in [50].
Using the sliding mesh for this complex application shows its potential to couple individual moving boundary-conforming meshes at common interfaces.

Single-Screw Extruder
In the following, we will apply the sliding mesh approach to a different application -still in extruding -but this time we consider single screw extruders (SSE). Similar to twin-screw extruders, single-screw extruders are built up by combining different screw geometries that each serve a different purpose. A selection of potential screw parts is shown in Fig. 13. Conveying elements (C) are used in order to transport the melt forward, whereas Maddock elements (B) serve as barriers that disperse the melt and decrease the pressure.
In contrast to twin-screw extruders, the conveying elements do not provide enough mixing. Thus, extra distributive mixing elements (E) have to be used. Finally, there is a section (G) that leads to the extrusion die, and other transition sections. An open research question is how to assemble the individual screw shapes in an optimal order, considering for example the mixing behavior or pressure loss, to name two potential design objectives [51]. The screw designs again feature very small gaps between screw and barrel, e.g., with a barrel diameter of 60 mm the smallest gap size between the conveying screw element and the barrel is only 0.3 mm. Thus, it is again of utmost importance to have a good discretization in those regions in order to capture all flow effects. This is where the sliding mesh approach comes into play. It allows to create individual boundary-conforming meshes for each individual screw section and to couple them at the common interface. Thus, it is not necessary to re-mesh in case one assembles the screws in a different way. In the following, we aim to demonstrate the usability of the sliding mesh approach to compute the flow field for four different screw combinations. The screw configurations are given in Table 6. The barrel diameter is 60 mm. As fluid, we consider corn syrup that can be modeled as a Newtonian fluid with density ρ = 1400kg/m 3 and viscosity η = 4.7Pa s. The screw rotates with rotational speed of ω s = 60 rpm. We set natural boundary conditions at the inlet and outlet. Thus, the flow rate will be mainly determined by the transport capacity of the conveying elements 1   2  3  4  combinations A-C1-B-D-E-F-G A-C1-C2-B-D-E-F-G A-B-C1-C2-D-E-F-G A-B-C1-D-E-F-G  elements  6,847,102  10,192,464  10,192,464  6,847,102   Table 6: Screw combinations based on Fig. 13 for four different screw combinations as well as the total number of elements for each resulting mesh.
in combination with the Maddock elements. For simplicity, we only simulate a steady-state flow without considering temperature effects.  The resulting mass flow rates are given in Table 7. Comparing configuration 1 and 4, as well as configuration 2 and 3, we can observe that there is a difference in the mass flow depending on the position of the Maddock element and the conveying elements. Adding an additional conveying element nearly doubles the mass flow rate, which is physically expected since the flow rate is mainly driven by the conveying elements. Furthermore, a factor that is slightly smaller than two makes sense since the mixing element also features a certain conveying behavior. The absolute increase of the mass flow rate by adding an additional conveying element is also very similar for configuration 2 (33.42 kg/h) and configuration 3 (32.28 kg/h). Fig. 14 shows the flow for all four configuration. The pressure field is shown on the screw surface and the velocity in z-direction on the plane xz-plane through the screw axis. The velocity as well as pressure contours are smooth across the non-matching interfaces. Furthermore, the mixing effects inside the mixing elements are clearly discernable by high negative velocities in z-direction. Additionally, the velocity magnitude reflects the different flow rates in case of one or two conveying elements. The pressure contours also demonstrate the pressure decrease over the Maddock element, which is mainly compensated by the conveying elements.

Conclusion
Within this work, we presented a stabilized finite element formulation for fluid and temperature equations that allows to couple sliding domains with boundary-conforming meshes at common interfaces. This approach allows to use highly-optimized boundary-conforming meshes for individual subdomains that can then be coupled weakly at their -possible only partially overlapping -common interface. We solve the non-isothermal, incompressible Navier-Stokes equations on the individual subdomains for Newtonian as well as for generalized Newtonain fluid models. The formulation is stabilized using residual based stabilization technique. The solution fields of each individual subdomain are coupled weakly. The interface coupling conditions are enforced using Nitsche's method. Additionally, the method imposes Dirichlet conditions weakly on boundaries that partly slide over rigid segments. We verified the method for the 2D Taylor-Green flow for viscous as well as convective flow regimes. Optimal convergence rates were obtained in a spatial convergence study, which suggests the correct scaling of our Nitsche stabilization terms. In a second validation step, a 3D test case considering the flow of plastic melt inside a twin-screw extruder kneading element with two discs was used to validate the method with respect to: (1) 3D cases, (2) generalized Newtonian models, and (3) the weak imposition of Dirichlet boundary conditions on the sliding interfaces. Again, optimal convergence rates for the interface coupling could be obtained in a spatial convergence study. First steps towards relevant industrial applications have been also made. We computed the time-dependent temperature and flow field of plastic melt inside a twin-screw extruder with several kneading and conveying elements. The SRMUM was used as mesh update method that allows to use boundary-conforming meshes for the individual rotating screw elements. The boundary conforming meshes were only coupled at the non-matching interface between individual kneading discs. As a second industrially relevant example, we computed the flow inside a single-screw extruder with different screw sections. Individual screw sections were meshed independently. We assembled the sections in four different ways. The resulting flow fields showed the expected behavior. This test case demonstrated the applicability of the presented approach to discrete optimization of single-screw extruders without any need for re-meshing. We believe that the presented method makes use of the benefits of both the boundary-conforming and unfitted approaches in a way tailored to the considered application cases. In the future, the method will be applied to more complex application cases, including realistic design optimization of extruders.