Analysis and optimization of mixing inside twin-screw extruders

Analysis and optimization of mixing inside twin-screw

extruders

Citation for published version (APA):

Sarhangi Fard, A. (2010). Analysis and optimization of mixing inside twin-screw extruders. Technische Universiteit Eindhoven. https://doi.org/10.6100/IR685851

DOI:

10.6100/IR685851

Document status and date: Published: 01/01/2010

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

Technische Universiteit Eindhoven, 2010.

A catalogue record is available from the Eindhoven University of Technology Library. ISBN: 978-90-386-2324-5

Subject heading: twin-screw extruders, mixing simulation, extended finite element, map-ping method

This thesis was prepared with the LA_{TEX 2ε documentation system.} Reproduction: University Press Facilities, Eindhoven, The Netherlands. Cover design: Arash Sarhangi Fard.

This research was conducted within the PEPTFlow project supported by the European Commission.

twin-screw extruders

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Eindhoven, op gezag van de

rector magnificus, prof.dr.ir. C.J. van Duijn, voor een commissie aangewezen door het College voor

Promoties in het openbaar te verdedigen op woensdag 29 september 2010 om 16. 00 uur

door

Arash Sarhangi Fard

prof.dr.ir. H.E.H. Meijer

Copromotoren: dr.ir. P.D. Anderson en

Summary ix

1 Introduction 1

2 Simulation of the viscous flow inside complex geometries 5

2.1 Introduction . . . 6

2.2 Mathematical formulation . . . 8

2.3 Numerical methods . . . 9

2.3.1 Lagrangian multiplier . . . 10

2.3.2 Non-conformal local mesh refinement . . . 12

2.3.3 The extended finite element method . . . 15

2.4 Results . . . 18 2.4.1 Geometry . . . 18 2.4.2 Boundary-fitted results . . . 20 2.4.3 FDM results . . . 20 2.4.4 XFEM results . . . 24 2.5 Conclusions . . . 26

3 3D simulation of the Stokes flow inside twin-screw extruders using XFEM 31 3.1 Introduction . . . 32

3.2 Mathematical formulation . . . 33

3.3 Numerical methods . . . 33

3.4 Results . . . 36

3.4.1 Validation . . . 36

3.4.2 Flow through the kneading elements of a twin-screw extruder . . . . 40

3.5 Conclusion . . . 48

4 Simulation of distributive mixing in twin-screw extruders using the mapping method 49 4.1 Introduction . . . 50

4.2 Basics of the mapping method . . . 53

4.2.1 The mapping matrix . . . 53

4.2.2 Concentration mapping . . . 55

4.2.3 Mixing measure . . . 55

4.3 Adaptation of the mapping method for TSE’s . . . 56

4.4 Application of the mapping method in TSE’s . . . 59

4.4.1 Velocity field . . . 59

4.4.2 Particle tracking . . . 62

4.4.3 Initial concentration distribution . . . 64

4.4.4 Screw assembly . . . 65

4.5 Conclusion . . . 66

5 Mixing and transport in conveying elements 67 5.1 Geometry . . . 67

5.2 Mixing analyses . . . 68

5.3 Influence of design parameters . . . 74

5.3.1 Gap size . . . 74

5.3.2 Pitch length . . . 76

5.4 Conclusion . . . 85

6 Mixing in kneading elements 87 6.1 Introduction . . . 87 6.2 Kneading elements . . . 88 6.2.1 Geometry . . . 88 6.2.2 Mixing . . . 90 6.3 SME elements . . . 97 6.3.1 Geometry . . . 98 6.3.2 Mixing . . . 98 6.4 Conclusion . . . 102 6.5 Reflection . . . 102

7 Conclusions and prospective 105

Bibliography 109

Samenvatting 115

Acknowledgments 117

The most commonly used mixing device in polymer processing is the closely inter- mesh-ing, co-rotating twin-screw extruder. The main goal of this thesis is mixing simulation inside twin-screw extruders using the mapping method. The thesis consists of two parts. First the computation of accurate three dimensional velocity fields using a combination of the finite element method, FEM, and the extended finite element method, XFEM, and second performing mixing analyses using particle tracking to define the mapping matri-ces and optimization of different configurations of screw elements by using the mapping method. The accuracy of mixing analyses depends on the accuracy of the velocity field. Determining the velocity field in twin-screw extruders, operating under realistic conditions, is still very challenging despite all progress made in simulations during last years. Two main difficulties arise when solving the balance equations. The first is that we have to deal with moving geometries, the screws, in a fixed barrel while the gap widths are extremely small and also change in time during rotation. Obtaining velocity fields and their deriva-tives, close to the moving screws is not trivial, but essential for further particle tracking analyses. The second challenge is to deal with a viscosity that is a function of position, shear rate and temperature, and can change orders of magnitude. In this work we present an extended finite element method to model non-Newtonian Stokes flow inside twin-screw extruders and demonstrate its accuracy and efficiency by systematically refining the mesh and compare with boundary-fitted results. Two-dimensional cross-sectional results for ve-locity and pressure are compared with full three-dimensional simulations. Residence time distributions are obtained for various screw designs. The analysis, and even better the optimization, of mixing based on particle tracking is far from trivial since no volumetric data are available. In this thesis we continue to use the mapping method, which does provide volumetric quantities, to quantitatively compare different screw layouts and find optimal designs. The mapping method has proven its merits for static mixers where, in previous work, all popular static mixers are evaluated and guidelines have resulted in the design of new mixers based on either compactness, minimum pressure drop, or accuracy in layer distribution. In this thesis we used to same approach to analyze mixing and com-pared several screw configurations and different types of screws. In particular, conveying elements are compared with different pitch length and gap width, while for kneading el-ements a variety of staggering angles have been selected and put next to each other. The flux-weighted intensity of segregation is used as a mixing measure to evaluate the different screw designs.

Introduction

Mixing is important in many polymer processing operations and significantly affects the properties of the final product and, moreover, its processability, energy consumption, and costs. In the design of new materials, blending and compounding techniques, helps to achieve the desired physical and mechanical properties, while blending is more cost-efficient than synthesizing new polymers. Polymer blends can show remarkable proper-ties directly related to their morphology, which results from phase deformation, breakup, coalescence, and retraction during processing. Therefore mixing of polymers is crucial. It determines the final morphology which accordingly strongly affects the properties of the blends.

Generally, the term mixing refers to operations that have a tendency to reduce non-uniformities or gradients; such as local differences in concentration, temperature, size of a dispersed phase, or other properties of materials. On a molecular level diffusion controls mixing, whereas on larger scales mixing is dominated by convection, the flow of materials [1]. Given the high viscosity of polymers, distributive mixing in laminar flows is the most relevant mechanism. Considering the lack of turbulence in these systems, the presence of chaotic advection determines the quality of distributive mixing, exponentially increasing the interfacial area between the different phases present in the blend.

Distributive mixing in polymer processing is performed in batch mixers or using continuous extruders. Among various types of mixing extruders, the co-rotating twin-screw extruder is the most widely used, given its flexibility. Co-rotating self-wiping twin-screw extruders are multi-purpose devices, and used for melting, pumping, reaction, degassing, filling and blending of the materials. For that, the screws are build up in a modular way applying many different types of screw elements. Conveying elements with different pitch angle, reverse or forward transporting, kneading elements with different staggering angles, also with positive, neutral or negative overall pitch, and other type of mixing elements are chosen with different sequence. Every specific mixing operation requires a dedicated screw design to perform and satisfy a particular mixing task. Given the importance to understand the evolution of mixing in the different screw configurations, simulation tools

are useful reducing the amount of trial and error in industrial practice.

In the past decades various investigations on simulating mixing in twin-screw extruders have been performed. Most of these studies are based on relatively straightforward anal-yses of the flow field (measuring e.g. the ratio of shear to elongational flows), on the computation of the residence time distribution and on the determination of particle trajec-tories [2]. In all these methods the tracked particles are mass-less, they have no volume either and, since mixing in well-designed mixers typically is of a chaotic nature, they are therefore of limited use only if one wants to investigate concentration distributions. The aim of this thesis is to study the dynamics of the mixing process in co-rotating twin-screw extruders using different types of screw elements and sequences. We aim at achieving a comparison of mixing efficiency changing geometrical parameters such as the gap size, the pitch length, the staggering angle, and/or introducing complete new mixer designs. This project consists of two parts: first obtaining an accurate velocity field, and second using this velocity field to analyse and optimize mixing.

The accuracy of any mixing analysis strongly depends on the accuracy of the velocity field. Three main difficulties arise when solving the balance equations. The first is that we have to deal with moving geometries, the screws, in a fixed barrel, and the second is that gap widths are extremely small and change in time during rotation. Obtaining a velocity field, and its derivative, very close to these tightly fitted moving screws is not trivial, but essential for further particle tracking analyses. A further challenge, the third problem, is to deal with a viscosity that is a function of position, shear rate and temperature, and that can change orders of magnitude in the high shear regions. Determining the velocity field in twin-screw extruders operating under realistic conditions is still very challenging, despite all progress made in simulation accuracy obtained during last couple of years. Ishikawa et al. [3] and Bravo et al. [4] used a classical finite element technique, re-meshing the geometry at each time step. Re-meshing of complex geometries, such as those in twin-screw extruders, is a tedious and time consuming process. It becomes worse, or nearly impossible, in case of structured meshes e.g. hexahedral elements in 3D meshes. Avalosse [5] presented a mesh superposition technique (MST) to address this problem. This technique is simple to implement but requires fine meshes in order to obtain a good description of the moving parts and to better conserve mass, especially in the case of geometries with small gaps. Recently, a distributed Lagrangian Multiplier/Fictitious Domain (FDM) technique was de-veloped with many similarities to MST method [6] which uses one, fixed, mesh. In all FDM based methods there is continuity of field variables across the interface between fluid and internal rigid body, since the fixed grid covers both the fluid and the rigid body domain, and there is no guarantee that the boundary of fluid elements aligns with the boundary of the rigid body. Therefore it is not possible to accurately predict discontinuities, or jumps, from fluid to solid along their interface for field variables such as velocity or pressure. This has a crucial effect on the accuracy of the computed velocity field in the gap regions where the pressure gradient is high and dominates the flow.

In this work we present an extended finite element method (XFEM) to model Stokes flow inside twin-screw extruders and we demonstrate its accuracy and efficiency by system-atically refining the mesh and compare results with those of boundary-fitted computa-tions which are feasible in a 2D geometry. We show that the accuracy of all FDM based methods has a fundamental problem, mainly in the gap regions, and even local mesh refinement does not help to increase its accuracy. In contrast XFEM predicts the flow in the gap regions close to the boundary-fitted results. In Chapter 2, the two-dimensional

cross-sectional results of velocity and pressure are presented comparing results of the boundary-fitted method with those of the refined fictitious domain method, and that of the XFEM. In Chapter 3 the full three-dimensional simulations of the velocity field are shown using the XFEM method for different screw elements.

Analysis, let it be optimization, of mixing based on particle tracking is far from trivial since no volumetric data are available. In this work we will continue on the application of the mapping method, which does provide volumetric quantities, to quantitatively compare different screw layouts and find optimal designs. The mapping method computes a distri-bution matrix that maps the color concentration distridistri-bution for a specific time, or specific length, to characterize mixing in a quantitative way. The calculation time for this method is efficient, since the mapping matrix is calculated only for one periodic unit and also for one periodic time [7].

The mapping method has proven its merits for static mixers where, in previous work, all popular static mixers are evaluated and guidelines have resulted to design new mixers based on either compactness of minimum pressure drop [8]. In this work we will use the same approach to analyze mixing and we will compare several screw configurations and different types of screws. The periodic property of the mapping method helps in handling mixing optimization for different screw elements, but also for many screw revolutions, since a twin-screw extruder is a time and space periodic device. Three-dimensional, volume mapping is done; therefore we capture mixing information over the different cross-sections of the extruder.

In particular, conveying elements are compared with different pitch length and gap width, while for kneading elements a variety of staggering angles have been put side by side. We also try to evaluate the mixing performance for the new designs of the screw elements such as the SME. The flux-weighted intensity of segregation is used a mixing measure to evaluate the different screw designs. Finally the performance of the different designs is compared from mixing and energy saving point of view. The details of the mapping method as applied to twin-screw extruders is described in Chapter 4 and the results of the mixing optimization for the conveying elements and mixing elements are shown in Chapters 5 and 6, respectively.

The work has resulted in a number of publications:

• Sarhangi Fard, A, Hulsen, MA, Meijer, HEH, Famili, NMH, and Anderson, PD. Adap-tive non-conformal mesh refinement and extended finite element method for viscous flow inside complex moving geometries. Int. J. Numer. Meth. Fl. 2010; submitted. • Sarhangi Fard, A, Hulsen, MA, Meijer, HEH, Famili, NMH, and Anderson, PD. Ex-tended Finite Element Method for Stokes flow inside complex 3D geometries with moving internal boundaries. Chem. Eng. Sci. 2010; in preparation.

• Sarhangi Fard, A, Meijer, HEH, and Anderson, PD. Simulation of distributive mixing in twin-screw extruders using the mapping method. Part 1: Transport elements.

Macromol. Theory Simul. 2010; in preparation.

• Sarhangi Fard, A, Meijer, HEH, and Anderson, PD. Simulation of distributive mix-ing in twin-screw extruders usmix-ing the mappmix-ing method. Part 2: Mixmix-ing elements.

Moreover, the author contributed to some papers that are not part of this thesis:

• Sarhangi Fard, A, Famili, NMH, and Anderson, PD. A new adaptation of mapping method to study mixing of multiphase flows in mixers with complex geometries.

Comput. Chem. Eng. 2008; 32(7):1471-1481.

• Diemert, J, Chilles, C, Colbert, J, Miri, T, Ingram, A, Sarhangi Fard, A, Anderson, PD. Characterization of residence time distributions in co-rotating twin-screw extrud-ers, using positron emission particle tracking experiments and flow simulation. Int.

Simulation of the viscous flow

inside complex geometries

1

Abstract

Two different techniques to analyze a non-Newtonian viscous flow in complex geome-tries with internal moving parts and narrow gaps are considered and compared. The first technique is a non-conforming mesh refinement technique based on the fictitious do-main method (FDM), and the second one is the extended finite element method (XFEM). The refinement technique uses one fixed reference mesh. To impose continuity across non-conforming regions, constraints using Lagrangian multipliers are used. The size of elements, locally in the high shear rate regions, is reduced to increase accuracy. FDM has limitations, therefore XFEM is applied to decouple the fluid from the internal moving rigid bodies. In XFEM the discontinuous field variables are captured using virtual degrees of freedom as enrichment and by applying special integration over the intersected ele-ments. The accuracy of the two methods is demonstrated by comparison with results of a boundary-fitted mesh applied to a 2D cross section of a twin-screw extruder. Compared to non-conforming FDM, XFEM shows a considerable improvement in accuracy around the rigid body, especially in the narrow gap regions.

1_{Reproduced from: A. S. Fard, M. A. Hulsen, H. E. H. Meijer, N. M. H. Famili and P. D. Anderson,(2010).} Adaptive non-conformal mesh refinement and extended finite element method for viscous flow inside complex moving geometries. Int. J. Numer. Meth. Fl., submitted.

2.1

Introduction

Numerical simulation of flow in complex geometries with internal moving rigid bodies, has shown great progress in the past decades. Modeling of moving discontinuities with the classical finite element method (or finite volume method) is cumbersome due to the need to update the mesh to match the geometry of moving rigid bodies. Re-meshing of com-plex geometries, such as a twin-screw extruder (TSE), is a tedious and time-consuming process. It becomes worse, or nearly impossible, in case of structured meshes e.g. hex-ahedral elements in 3D meshes. To overcome these problems, several numerical tech-niques were developed. We can classify these methods in two categories: moving mesh and fixed mesh.

One of the important moving-mesh methods is the Arbitrary Lagrangian Eulerian (ALE) technique, which is frequently used to simulate the interaction of flows with a flexible structural domain [9, 10]. This method is based on using a deforming grid for the fluid. The grid deforms with the flexible structure at the interface and then the grid deformation is extended into the fluid field. Since the boundary of the fluid mesh completely aligns with the boundary of the rigid body, the quality of the deformed mesh is preserved, and the accuracy of the ALE method is higher than any fixed mesh method. ALE methods have some limitations if the internal moving part has a rigid complex geometry, since the fluid mesh highly deforms and becomes distorted. Sometimes re-meshing helps to overcome this problem but in case of 3D problems it is really complicated, or impossible, to deform the finite element mesh with flow advection (e.g. the 3D geometry of a TSE). On the other hand re-meshing of a deformed mesh is, from a computational time view, expensive. To circumvent deforming mesh problems, fixed mesh methods were introduced based on a reference mesh and immersing boundaries of moving parts in the fixed mesh. The original immersed boundary (IB) method was introduced by Peskin [11] and uses a fixed mesh for the computational domain, that is intersected by the fluid/solid interface. Using a level-set function, the additional interface nodes are added to the computational domain with additional degrees of freedom. A special treatment is then applied to the solid region and the interface degrees of freedom [12]. Avalosse [5] presented a mesh superposition technique (MST) applied to a twin-screw extruder. With MST the whole computational do-main is covered with one fixed mesh without considering the internal moving rigid bodies and one separate dynamic mesh for the moving parts. At each time the position of the moving part mesh is updated. With superposition of these two meshes, the kinematics of the dynamic mesh are imposed on the static mesh using penalty constraints. This method is easy to implement but to avoid diffusion errors, precise mesh refinement at the interface of fluid-rigid body is required [13].

Recently, a distributed Lagrangian Multiplier/Fictitious Domain (FDM) technique was de-veloped with many similarities to MST or IB methods [6, 14]. With FDM, a fixed mesh for the computational domain, including internal moving parts, is created. The kinemat-ics of the rigid-body motion is enforced with a constraint on the conservation equations. This constraint is imposed via Lagrangian multipliers. The rigid-body boundary has been discretized through collocation points. The rigid-body motion is imposed on (evenly dis-tributed) points on the solid boundary (a line in 2D or a surface in 3D). D’Avino and Hulsen [15] uses weak constraints for the rigid-body boundary discretization. The par-ticle boundary is divided into elements where the rigid-body motion is imposed in a weak sense. In this work we will compare the collocation points approach with the weak

con-straint approach.

The main challenge of FDM is prescribing precisely the enforced velocity of the moving rigid body on the fluid flow. The accuracy becomes critical in high shear rate regions. For example, in a TSE there are narrow gaps (screw-screw and screw-barrel) where the shear rate is orders of magnitude higher than in other parts of the flow regions. The accuracy of FDM in high shear rate regions can be increased by locally decreasing the element size of the fixed mesh (mesh refinement technique) across the interface of rigid body-fluid. Bertrand et al. used standard and non-standard mesh refinement techniques for the quadratic Crouzeix-Raviart triangular element,P2+−P1for a 2D twin-screw problem

based on the FDM [16], and later extended their method to 3D tetrahedron elements [17]. The challenge of local mesh-refinement is to impose continuity between refined elements and coarse elements. Bertrand et al. modified the shape functions associated with the non-conforming nodes. This method is easy to implement as it only requires a few mod-ifications of the shape functions and the assembling procedure. Another refinement ap-proach is using mortar element techniques. This method is introduced by Bernardi et al. [18, 19]. The application of this method is also developed for spectral elements [20]. Mortar element methods solve the coupling of the new refined elements with coarse ele-ments via adding extra constraint equations to the stiffness matrix. Belgacem [21] used a Lagrangian multiplier for coupling the non-conforming regions with the rest of the com-putational domain. This method preserves the quality of the reference mesh. Its main drawback is the increasing number of constraints that may lead to an ill-conditioned stiff-ness matrix.

In all FDM based methods there is continuity of field variables across the interface be-tween fluid and internal rigid body, since the fixed grid covers both the fluid and the rigid body domain and there is no guarantee that the boundary of fluid elements aligns with the boundary of the rigid body. Therefore it is not possible to accurately predict disconti-nuities or jumps from fluid to solid along their interface for field variables such as velocity or pressure. Capturing of discontinuous field variables such as pressure, has a crucial effect on accuracy in the gap region where the pressure gradient is high and dominates the flow.

The Extended Finite Element Method (XFEM) is originally introduced for crack prob-lems to predict the discontinuity of field variables such as displacements and stress in solids [22, 23]. Recently XFEM was further developed to include applications in fluid me-chanics. Wagner et al. [24] used XFEM to simulate problems of rigid particles in Stokes flows using analytic solutions as the partition of unity enrichment. Gerstenberger and Wall [25] applied XFEM to problems of fluid-structure interaction in Newtonian fluids us-ing a Heaviside function as the enrichment. Choi et al. interpreted the XFEM enrich-ment scheme using virtual degrees of freedoms and applied the method to particulate viscoelastic flows [26]. Chessa and Belytschko [27] applied XFEM to two-phase fluids. In this work we present a new strategy for a non-conforming mesh refinement method based on FDM for rectangular Taylor-Hood elements and its adaptation for complex ge-ometries with internal moving parts and narrow gap regions. The local mesh refinement technique presented is based on one fixed reference mesh. At each time step the size of the intersected elements containing the boundary of the moving parts is decreased by splitting elements. To ensure continuity of non-conforming refined elements with coarse elements we use extra constraints, implemented by using a Lagrangian multiplier. The quality of the refined mesh is the same as that in the reference mesh. A 2D cross

sec-Ω d Γ 2D cross-section s Γ

Figure 2.1: 2D cross section of a self-wiping co-rotating twin-screw extruder with narrow gaps between screw-screw and screw-barrel.

tion of a twin-screw extruder with the typical narrow gap regions between screws and the barrel, and between the screws themselves, is chosen as a test case. The flow inside this geometry is the result of a combined drag and pressure flow, with high shear rates in the gap regions. Next, we employ the XFEM enrichment scheme using virtual degrees of freedom. We compare the results of the non-conforming mesh refinement using FDM with XFEM to explore the effect of discontinuities of physical variables.

2.2

Mathematical formulation

The polymer melt inside the flow domainΩ with boundariesΓ (Figure 2.1) is modelled as an isothermal and incompressible fluid. The Reynolds number is low and inertia will be neglected. Also viscoelastic effects are ignored and the fluid is described using a generalized Newtonian constitutive equation. The flow is governed by the momentum and the continuity equations:

−∇ ⋅_{τ + ∇p}= 0 _{inΩ, (2.1)}

∇ ⋅_v= 0 _{inΩ, (2.2)}

where v is the velocity vector andpis the pressure. The relation between the extra stress tensorτ and the symmetric part of the velocity gradientDis given by:

τ = 2η(IID)D, (2.3)

whereD= 1₂[∇v + (∇v)T]is the rate of deformation tensor andIIDis the second invariant ofD. The Carreau model is used:

η( ˙γ) = η∞+(η0−η∞)(1 + (κ ˙γ)2)(n−1)/2, (2.4)

whereη_∞,η0,κand nare parameters of the Carreau model [28] and the shear rateγ˙ is

related to the second invariant ofD:

˙ γ=

√ 1

Dirichlet boundary conditions for the velocity are imposed onΓ :

v= vd onΓd, (2.6)

v= vs onΓs, (2.7)

wherevd andvs are the velocities on moving (Γd) and stationary (Γs) boundaries,

respec-tively.

2.3

Numerical methods

The governing equations are partial differential equations that can be solved with the finite element discretisation method [29]. The weak form of the conservation equations is obtained by multiplying Equations 4.14 and 2.2 with test functionsν and qrespectively, and partially integrate the momentum balance:

((∇ν)T_{, τ) − (∇ ⋅ν, p) = 0 ∀ν ∈ (H}1

0(Ω))2, (2.8)

(q, ∇ ⋅ v) = 0 ∀q ∈ L2_{(Ω), (2.9)}

where(f, g) = ∫_Ωf ∶ g dA.

The computational domain is discretized using a finite element mesh. The field variables,

vandp, are discretized as:

vh(x) = nv ∑ i=1 φ_i(x)v_i ∀_x∈ Ω, (2.10) ph(x) = np ∑ i=1 ψ_i(x)p_i ∀_x∈ Ω, (2.11)

wherevh and ph are the interpolated functions of the velocity and pressure in Ω, nv is

the number of velocity nodal points, np is the number of pressure nodal points, φi(x)

and ψ_i(x) are the shape functions of velocity and pressure inΩ, and vi and pi are the

nodal values of velocity and pressure, respectively. Rectangular Taylor-Hood elements are used; nine nodes (corners, mid-sides and center) for velocity with bi-quadratic shape functions and four nodes (corners) for pressure with bi-linear shape functions. Since the viscosity is shear-rate dependent, the non-linear part of the equations is linearized using a Picard iteration. The numerical procedure is described in 2D; generalization to 3D is straightforward.

In the case of geometries with moving internal parts, like for example twin-screw extruders (Figures 2.2a-2.2c), the flow domain changes. To avoid the need of re-generation of the mesh for each position of the moving part, the fictitious domain method (FDM) is applied [6, 14]. In FDM, the kinematic condition of the moving parts is enforced using constraints that are implemented with a Lagrangian multiplier (Figure 2.2d):

((∇ν)T_{, τ) − (∇ ⋅ ν, p) + (ν, λ)}

Γ_d = 0, (2.12)

(q, ∇ ⋅ v) = 0, (2.13)

(a) (b)

(c) (d)

Figure 2.2: Boundary fitted finite element mesh for 2D cross section of twin-screw ex-truder at three screws positions: (a) 0○, (b) 49○, (c)90○ and (d) reference mesh for fictitious domain method with constraint mesh of screws for0○ (○) and40.5○(◁) screw positions.

whereµ is a test function for the Lagrangian multiplier λ and (., .)Γ_d denotes the proper

inner products on the interface of the rigid body.

2.3.1 Lagrangian multiplier

Two approaches exist to impose the kinematics of a moving object: the first is using collocation points (control points) to depict the geometry of the moving part (Figure 2.3b). In this approach the constraint is imposed for each collocation point, therefore Equation (2.14) becomes:

µ_k⋅(vk−_vk_d) = 0, _k= 1, ⋯, J, ∀µ_k, (2.15) where J is the number of collocation points, µ_k is a arbitrary value (test value) at collo-cation pointk,vk is the interpolated velocity at the coordinate of collocation pointk, and

vk_d is the velocity imposed on the moving object at the collocation point k. The number of collocation points should be chosen such that it allows to cover all elements that are intersected with the moving object. To avoid an over-constrained system, the maximum number of collocation points per element should not be more than two [17]. The ap-proach is easy to implement, but when applying local mesh-refinement, J should also dynamically change, dependent on new intersected elements. The approach enforces constraints only in discrete points (collocation points) and there is no continuity between values of imposed constraints along the object-fluid interface. This becomes more sig-nificant if a higher-order integration scheme is used (e.g. XFEM). The other issue is that there is no guarantee to have collocation points for those elements which have a very small overlap with the rigid body.

To overcome these problems we used an alternative approach: a weak constraint on the object-fluid interface, instead of on the individual collocation points (Figure 2.3c). Based

W

H

V

Figure 2.3: (a) Schematic of a lid-driven cavity (H= W) with a fixed object in the center

(R/W = 1/8) and moving upper wall. Typical mesh (in the computations, a

much finer mesh is used); the non-rotating object is considered as a constraint on the boundaryΓ_d: (b) with collocation points; (c) with weak constraint.

on Equation 2.14 the constraints are enforced in a weak way by integration over the entire length of the interface mesh [15]:

∫_Γ

d

µ ⋅(v − v_d) ds = 0 ∀µ ∈ L2(Γ_d). (2.16) With this method the object-fluid interface is described with a mesh that possesses a topology. We use a linear shape function to interpolate the Lagrangian multiplier over the discretized interface: λh(x) = nλ ∑ i=1 Qi(x)λi. (2.17)

Herex represents the coordinates of the point on the boundary of the rigid body (Γ_d) on which we wish to find an approximate (interpolated) value of the Lagrangian multiplier

λ_h(x),i is the node index, nλ is the total number of nodes of the interface mesh, λi are

the nodal values of the Lagrangian multiplier, andQi(x)is the shape function associated

with nodeiinΓd.

In the weak constraint approach, the constraints are continuously distributed over the interface of rigid body. Therefore, a guarantee exists that with a sufficient number of integration points, the intersected elements are influenced by a constraint of the moving object. Similar to the collocation points approach, the number of elements of the interface mesh should be balanced with the number of intersected elements of the finite element mesh, but it is much less critical here.

To show the performance of the weak constraint approach, and to compare it with the collocation points approach, we first use the simple lid-driven cavity flow as a test case (Figure 2.3a). The object (a disk) in the center of cavity is fixed and the upper boundary of the cavity is dragged with dimensionless speedVd (Vd= 1). The cavity mesh consists

of 6400 (80 × 80) quadrilateral elements. The presence of the object is described by using Lagrangian multipliers and the problem is solved with collocation points and weak constraints. In the weak constraint method, the boundary of the disk is divided into one-dimensional quadrilateral elements where the constraints are imposed in an integral way. The integration is performed by using subintervals through a mid-point integration rule and results are obtained for a different number of integration points per element (1, 10

0 45 90 135 180 225 270 315 360 −0.005 0 0.005 0.01 0.015 0.02 0.025

θ (degree)

collocation points constraint

weak constraint with 30 integration points weak constraint with 10 integration points weak constraint with 1 integration point

V

/V

d

θ

Figure 2.4: Comparison of computed non-dimensional tangential velocities (V/Vd) at some sampling points along the fluid-object interface, using the collocation point approach and the weak constraint approach with different number of intergration points.

and 30 points), see Figure 2.4. It is observed that for the collocation point constraint approach the velocity at the interface for the positions of the collocation points is zero but between collocation points, for any sampling position, it is non-zero and varies. For the weak constraint approach this error is much lower and decreases to zero when the number of integration points is sufficiently high.

2.3.2 Non-conformal local mesh refinement

Moving internal boundaries and narrow gaps are two important issues in the simulation of the velocity field inside complicated geometries. In the previous section it was described that to overcome regeneration of the mesh for geometries with internal moving parts, the velocity at the boundary of the moving parts can be enforced using a Lagrangian multiplier

λ. To implement this method, the boundaries of the moving parts are described with a surface mesh, and the accuracy of the fictitious domain method depends on the number of elements intersected with moving boundaries and the degree of approximation of λ. The resolution depends on the number of elements used in the moving boundary mesh (or the number of collocation points), but with an increasing number of constraints, the number of intersected elements of the main domain should also be increased. In geometries with narrow gaps, high velocity gradients occur and their positions change in time. To obtain adequate accuracy in these regions, the size of the elements there should be small. Local mesh refinement is efficient but it leads to a non-conformal mesh. Here we use the mortar element method [21] and an automatic mesh refinement process is designed which relies on one fixed reference mesh. First criteria have to be defined to detect target elements of the reference mesh for refinement. In the fictitious domain method an algorithm is

e Ω 1 e Γ 2 e Γ 3 e Γ 4 e Γ (a) 1 e Ω e2 Ω 3 e Ω e4 Ω 1 e Γ 2 e Γ 3 e Γ 4 e Γ (b) 1 e Γ 2 e Γ 3 e Γ 4 e Γ (c)

Figure 2.5: Non-conformal subdivision of a 2D quadrilateral element: (a) target element

Ωe; (b)2× 2subdivision; (c)3× 3subdivision. Coupled nodes are shown by∎ and discontinuous nodes by●.

used to find all elements that intersect the internal moving boundaryΓ_d. Those elements are located in the narrow gap regions, thus they are refinement targets. Next we subdi-vide the target elements into small elements. For example in Figure 2.5 for quadrilateral nine-nodes elements, the splitting procedure is shown for 2 × 2 (Figure 2.5b) and 3 × 3

(Figure 2.5c) subdivisions. All sub-elements and target elements have the same order, although it is possible to generate lower or higher order sub-elements. The challenge of non-conformal refinement is to solve issues related to the discontinuity of the field vari-ables such as velocity or pressure between the refined elements and the coarse mesh, since the new nodes of sub-elements do not belong to coarse elements. For example if the elementeis divided into 4 sub-elements (refinement order 2), then:

Ωe= Ωe₁∪_Ωe₂∪_Ωe₃∪_Ωe₄. (2.18)

All non-coupled new nodes are located at the boundary of Ωe (Γ_ie,_i = {1, 2, 3, 4}) and these non-conforming boundaries should be coupled to the reference mesh. Since the refinement of elements is processed one by one, it is more efficient from a computa-tional point of view to connect all new sub-elements that have a same order of refinement (Figure 2.6a).

To ensure the continuity of field variables between refined elements with neighboring coarse elements, we use a weak constraint implemented with a Lagrangian multiplier over the non-conforming interface. In a general geometrically non-conforming domain (Figure 2.6b) the mortars areγ_i ∈ {Γ1,2, Γ1,3, Γ2,3}. We have to minimize the difference in

function values (e.g. velocity) across each mortar (γ_i) using a Lagrangian multiplier (λ) and a test functionµin a weak constraint:

∫_Γ 1,2 µ ⋅(v1−v2) ds = 0 (2.19) ∫_Γ 1,3 µ ⋅(v1−v3) ds = 0 (2.20) ∫_Γ 2,3 µ ⋅(v2−v3) ds = 0, (2.21)

and the corresponding (symmetric) contribution to the momentum balance. Discretization ofλis defined on the boundaries of the refined elements.

1 Γ 2 Γ 3 Γ 4 Γ 5 Γ 6 Γ 7 Γ (a) (b)

Figure 2.6: (a) Connecting of side-by-side sub-elements with the same order of refine-ment; (b) Mortar decomposition in a geometrically non-conforming case, three mortars. 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 (a) (b) (c) (d)

Figure 2.7: Local mesh refinement technique is applied for a cavity with rotating rigid body in the center. Accuracy of non-conforming mesh refinement technique is rep-resented with calculated average shear rate inside the rigid body; (a) no re-finement γ¯˙ = 0.66 s−1, (b) 2× 2 subdividing γ¯˙ = 0.15 s−1, (c) 3× 3 - 2× 2

×

×

×

×

×

×

×

_×

Ω

Ω

e

e

e

e

Figure 2.8: (a) One side discontinuity problem between fluid domain (inside)Ω+and rigid body domain (outside)Ω−. Non-intersected elements fully insideΩ+are fully integrated, intersected elements (gray color) are integrated on the inside part (Ω+) and the elements (×) fully outside (Ω−) are ignored. Virtual degrees of freedom in nodes outside but belonging to intersected elements are shown with●. (b) Subdivision of elements: no subdivision for fully inside elemente1, two level quad-tree subdivision for elemente2, three level quad-tree subdivi-sion for elemente3, and four level quad-tree subdivision for elemente4. Final quadrilateral sub-elements are decomposed to small triangles to describe the interface.

This mesh refinement strategy is first applied to simulate a fluid flow inside a cavity with stationary boundaries and one rotating circular object (speed 1 rev/s) in its center. The geometry and fluid properties are the same as in Figure 2.3 except for the radius of the object (R/W = 2.25/8). The reference mesh consists of 15 × 15 quadrilateral elements with nine nodes (Figure 2.7a). Different mesh refinements are applied. In the first one we subdivide those elements that intersect with the object to2 × 2sub-elements (Figure 2.7b). Next is subdividing intersected elements to3 × 3, and their neighboring elements to2 × 2

sub-elements (Figure 2.7c). The last one is subdividing intersected elements to 4 × 4

sub-elements, and neighboring elements to 3 × 3 and 2 × 2 sub-elements (Figure 2.7d). The computed shear rate is used to investigate the performance of these different mesh refinements see Figure 2.7. With increasing the refinement steps, the calculated mean shear-rate inside the object decreases to zero.

2.3.3 The extended finite element method

In the fictitious domain method, the fluid flow problem is solved with a finite element mesh which includes also the solid moving part. This means that the boundary of the solid body is not aligned with the fluid element boundaries; intersected elements exist which are partly inside the fluid and partly inside the rigid body. As a result, there is a non-physical continuity for field variables such as velocity and pressure across the interface of the fluid-rigid body along these intersected elements. To include jumps for field variables along the interface, we use the extended finite element method (XFEM). For a one-sided discontinuity problem (Figure 2.8) the intend is solving the fluid flow (Ω+) only.

The variation of each unknown continuous field variable, such as f, overΩcan be inter-polated using shape functions:

f(x) ≈

nf

∑ i=1

Ni(x)fi, (2.22)

wherexrepresents the coordinates of the point inΩon which we wish to find an approx-imate (interpolated) value for the function f, i is the node index, nf is the total number

of nodes in Ω, fi are the nodal values of the function f, and Ni(x) is the shape

func-tion associated with nodei in Ω. To impose the jump for field variables such as f, the approximation function of f inΩcan be redefined as:

f(x) ≈

nf

∑ i=1

Ni(x)H(s)fi, (2.23)

whereH(s)is the Heaviside function defined with a scalar level set functions: H(s) = { +1 if s≥ 0 0 if s< 0, (2.24) s(x)⎧⎪⎪⎪⎨⎪⎪⎪ ⎩ > 0 if x∈ Ω+ < 0 if x∈ Ω− = 0 ifxis on the interface. (2.25)

The level set functions(x) is also used to determine whether elements are fully inside, fully outside or intersected with the fluid-rigid body interface. The definition of the level set function is based on the shape of the rigid-body boundary. For a regular geometries with boundaries that can be described with simple mathematical equations such as for a cylinder, it is straightforward to defines(x)analytically. For non-regular or complicated geometries,s(x)is defined approximately or numerically.

If we have the situation like in Figure 2.8a, elements (×) fully outside (Ω−) are ignored, elements fully inside (Ω+) are fully integrated as before, but elements (shaded elements) that are partially inside and outside need only to be integrated on the inside part see the Figure 2.8b. The degrees of freedom in the inside nodes (∎) have the same meaning as before, they are standard degrees of freedom. The nodes that are connected to fully outside elements are removed from the solution and fully outside elements are not as-sembled in the final solution matrix and field variables in their nodal points are set to zero as Dirichlet boundary condition. The degrees of freedom in nodes outside of (Ω+) but belonging to intersected elements are virtual degrees of freedom (●). They do not define the field variable at that point, but define the spatial distribution of the field variable in the part of the element that is inside the body. This is similar to virtual or ghost grid points in the finite volume method. The value in the ghost nodes is just considered for interpolating values at the fluid domain part of an intersected element. Using the XFEM method, the velocity and pressure are redefined as follows:

vh(x) = nv ∑ i=1 φ_i(x)H(s)v_i (2.26) ph(x) = np ∑ i=1 ψ_i(x)H(s)p_i (2.27)

+

Ω

Ω

a

a

Ω

Ω

ˆ

n

Figure 2.9: (a) An intersected element: a1is the area of the inside part of the element and

a2is the area of the outside part of the element; (b) Moving mesh (shown with

dashed line) in the direction normal to the screw surface.

An important part of XFEM is the numerical integration scheme. The intersected elements are subdivided into sub-domains and integration is performed over these sub-domains. The integration on sub-domains of an element determines the robustness and accuracy of XFEM. The XFEM integration scheme for quadrilateral elements is shown in Figure 2.8b. The first step is to create a quad-tree; using the level set function, and the intersected ele-ments are determined. Then the quad-intersected eleele-ments are subdivided to small quad elements (1-level subdividing), new intersected sub-elements are determined and further sub-dividing them to small sub-elements (2-level subdividing) follows. With increasing the level of quad-subdivision the integration area (or volume) is much better defined by the interface geometry which increases the accuracy of integration. It also increases the CPU time, so finding the optimum number for the subdividing level is important. In this work a 5-level subdivision is used. After creating the quad-tree, the smallest quad inter-sected sub-domains are subdivided into small triangles. The algorithm of triangulation using a level set function is well described in [30]. Now after domain subdividing, the flow part aligns with the interface and numerical integration can be done only on sub-domains (triangles or quadrilateral sub-sub-domains) that are fully inside. We used a3 × 3

Gauss integration rule for the quadrilateral elements and sub-elements, and a 6-points Gauss integration for triangles.

Also an internal Dirichlet boundary condition onΓ_d is imposed using a weak constraint implemented with a Lagrangian multiplier, similar to FDM.

In Figure 2.9a an intersected element is shown. If the integration area (a1) is very small

compared to the element area, the final solution matrix becomes ill-conditioned. This situation becomes worse when iterations on the solution are needed, e.g for non-linear problems and a divergence of the solution can occur. To avoid this problem, we used a method introduced by Choi et al. [26] that moves the mesh inside the domain slightly along a direction normal to the interface of the fluid-rigid body (Figure 2.9b), moving nodes in-side the fluid, while nodes inin-side the rigid body are fixed. By solving the Poisson equation,

the velocity of nodesun can be obtained:

∇2_un= 0 inΩ (2.28)

un = 0 onΓs (2.29)

un = αen onΓd (2.30)

wherenis the outwardly directed unit normal vector on the rigid body surface andαeis a coefficient to control the value ofnin elemente with a default value of unity. Each mesh pointxn moves according to the following advection equations:

dxn

dt = {

un if xn∈ Ω+

0 if xn∈ Ω− (2.31)

xn(tn= 0) = xn,0 (2.32)

With xn,0 the original position of the unmoved mesh nodes. A fourth-order Adams-Bashforth method is used to solve the system (2.31). The nodal points outside the rigid body are moved until each integration area of intersected elements is larger than 3% of the entire element area (a1/(a1+a2) > 0.03). Moving the mesh may distort some elements.

Displacing of the mesh is controlled by the time step∆tof Adams-Bashforth method, and value ofαe. To avoid large local distortions of elements the following procedure is applied:

• A node in-between element vertices is set to the midpoint of vertex nodes after each Adams-Bashforth time step.

• The coefficient αe for elements located in the gap regions is set to less than unity (0.5 < αe < 0.9) and in the other regions is set to unity. Sometimes if the gap is extremely narrow, a negative value (e.g. αe = −0.05) is used, moving the mesh in opposite direction.

• We use small values of∆t,∆t= 0.001.

2.4

Results

2.4.1 Geometry

To validate the numerical methods introduced, the cross section of a twin-screw extruder is selected as a test geometry (Figure 2.10). It consists of two moving internal parts. Nar-row gaps are present between screw-screw and screw-barrel. For a given screw radius (Rs), centerline distance (Cl), screw-screw gap distance (δs), barrel-screw gap distance

(δ_b) and number of parallel channels (e) it has a unique shape because of the requirement that one screw must, in any position, wipe the other one. A number of different descrip-tions for the screw geometry of co-rotating twin-screw extruders exist [31, 32, 33]. We use

x y l C s Rδb s δ α ( ) l θ * ( )

r

Figure 2.10: Geometric parameters of self-wiping twin-screw extruder.

a slightly adapted version of the description of a self-wiping profile provided by Booy [31]:

l(θ) = ⎧⎪⎪⎪⎪ ⎪⎪⎪⎪ ⎨⎪⎪⎪ ⎪⎪⎪⎪⎪ ⎩ Rs(ρc−1) 0≤ θ ≤ 1₂α Rs[ √ ρ2c−sin2(θ −α2) − cos(θ − α 2)] 1 2α≤ θ ≤ 1 2α + 2β Rs 12α + 2β≤ θ ≤ 1 eπ , (2.33)

wherel(θ)is the distance from the center to the edge of each corresponding screw. The geometrical relations of the parameters are:

ρc= Cl −_δs Rs , (2.34) β=arccos(ρc 2 ) , (2.35) α= π e −2β. (2.36)

The values of the parameters are listed in Table 2.1.

Based on Equation 2.33 a level set function s(x) is defined to determine whether the points are inside(s(x) < 0), outside(s(x) > 0)or on the screw surface(s(x) = 0):

s(x) = r(θ) = ∣OP∣ − l(θ) ∀P ∈ Ω, (2.37)

where r(θ) is a linear distance function of point P with coordinate x from the center of the screws (O or O´ in Figure 2.10). If the element intersects with one screw (e.g. left screw), in the XFEM integration scheme the related reference point for that screw (e.g. O) is considered. For those elements that intersect with both screws (frequently in the intermeshing region), elements are subdivided into small elements until their subdivided elements intersect with just one screw. Note,s(x)is not the real signed distance function for the screw surface and it defines a linear signed distance (shortest distance) from the reference point (OorO´).

We use a shear-thinning fluid in order to show the differences in numerical methods pro-posed, since in the intermeshing and gaps regions shear rates are high, and strongly influencing the viscosity. The values of rheological parameters for the Carreau model are listed in Table 2.2.

Table 2.1: Screw geometry parameters

Number of flights (e) 2

Screw radius (Rs) 15.275mm

Center distance of screws (Cl) 26.2mm

Clearance between screws (δs) 0.20mm

Clearance screw-barrel (δ_b) 0.15mm

Table 2.2: Rheological parameters for the Carreau model.

η0 1290 Pa⋅s

η∞ 0.0 Pa⋅s

n 0.559

-κ 0.112 s

2.4.2 Boundary-fitted results

To compare results, we used a classical finite element method with a single bound-ary fitted mesh that precisely takes into account the correct interface conditions. Since the geometry is changing with time, we just consider a few orientations. Fig-ure 2.11 shows typical boundary fitted (BF) meshes for selected orientations {θ =

0○, 22.5○, 90○, 112.5○and157.5○}, Figure 2.12 shows the details of this mesh. The ve-locity values are calculated forθ= 0○ and for screw rotation speedω= 1rev/s along the intermeshing position x = 10.8246mm (the origin of Cartesian coordinates is located in the center of left screw), using different number of nodes, see Figure 2.13. With decreas-ing element size, the values of the tangential velocity converge. A BF mesh with 118877 nodes was selected as a reference result to make the comparison between all numerical methods introduced.

2.4.3 FDM results

Non-conforming local mesh refinement was described in the previous section to increase the accuracy of the fictitious domain method by reducing the size of intersected elements. We apply non-conforming FDM on two reference meshes F1 and F5. For each reference mesh three refinement strategies are used: 2 × 2 subdividing for intersected element,

3 × 3subdividing for intersected elements plus2 × 2subdividing for their neighboring ele-ments and4 × 4subdividing for intersected elements plus2 × 2subdividing for their neigh-boring elements. Mesh information is summarized in Table 2.3. As a demonstration, a simple quadrilateral structured mesh, represented in Figure 2.14, is refined with the three schemes using the non-conforming mesh refinement technique and the shafts of the screws are considered hollow, to reduce the number of elements.

For one screw revolution, increments of4.5○ are chosen for the time steps between two consecutive screw configurations, corresponding to a time interval of 0.0125 seconds for a counter clockwise rotation speed of 1 rev/s. Since the geometry of the two-lobe screw elements is symmetric, only one half of the complete cycle is considered in the

(a)θ= 0○

+

(b)θ= 22.5○ (c)θ= 45○

(d)θ= 90○ (e)θ= 112.5○ (f)θ= 157.5○

Figure 2.11: Typical BF meshes with quadrilateral nine-node elements for different screw orientation.

Figure 2.12: Typical BF mesh of twin-screw extruder with98431nodes and mesh quality in the gap regions.

Table 2.3: Meshes used for simulations of FDM and mesh refinement.

nelem1 nnodes2 RS3 nelem14 nnodes15 S6 CPU time (s)

F1 29791 120921 - - - 0.138 123 F2 29791 120921 2 33019 140883 0.069 426 F3 29791 120921 3 − 2 45103 201285 0.049 1745 F4 29791 120921 4 − 2 53398 236815 0.034 2220 F5 52726 213257 - - - 0.103 238 F6 52726 213257 2 57478 241769 0.051 774 F7 52726 213257 3 − 2 73249 321075 0.034 4469 F8 52726 213257 4 − 2 84337 368595 0.026 5602

1_{Number of elements for reference mesh} 2_{Number of nodes for reference mesh} 3_{Refinement scheme}

4_{Number of element after refinement} 5_{Number of nodes after refinement} 6_{Size of intersected elements,}

S=(circumference of screw)/(number of intersected ele-ments), circumference of screw≈ 81.67mm

−0.80 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 5 10 15 20 25 BF1=62521 nodes BF2=98431 nodes BF =118677 nodes BF3=171077 nodes −5 −4 −3 −2 −1 0 1 2 3 4 5 x 10−3 24.75 24.8 24.85 24.9 y/Rs V /V p BF1 BF2 BF _BF3 θ=0 x = 10:8246 mm

Figure 2.13: Comparison of non-dimensional tangential velocity _VV

p (Vp= 2πRsω) along

the inter-meshing line (x= 10.8246mm) for different size of BF mesh. With increasing number of nodes the results converge to the same values.

(a) (b)θ= 99○

(c)θ= 310.5○ (d)θ= 54○

Figure 2.14: Example of a mesh F5 used for fictitious domain method and mesh refine-ment; (a)the reference mesh F5, (b)2× 2subdividing for intersected element F6, (c)3× 3subdividing for intersected elements and2× 2for their neighbor-ing elements F7, (d)4× 4subdividing for intersected elements and 2× 2for their neighboring elements F8.

−1 −0.95 −0.9 −0.85 −0.8 −0.75 −0.7 0 0.2 0.4 0.6 0.8 1

x/Rb

V

/V

p

Figure 2.15: The non-dimensional tangential velocity (_VV

p) as a function of

x

Rb over the line

y = 0 mm at the left side screw for the screw positionθ = 0○, using differ-ent meshes and differdiffer-ent refinemdiffer-ent scheme and comparison with boundary fitted result (BF).

simulations.

Figure 2.15 shows the non-dimensional tangential velocity (_VV

p) as function of

x

Rb over the line y= 0, x < 0for screw configurationθ= 0○. Results are shown for meshes in Table 2.3 and compared with the boundary fitted result. Generally with increasing total number of elements, the velocity values move somewhat toward the BF result, but the error still remains relatively large (with 6% average error at best for F8). To explore the performance of the local refinement technique we also show CPU time and size of intersected elements S, in Table 2.3. With increasing the subdividing, the size of intersected elements locally decreases, and the accuracy of the results improves. For example the accuracy of finest non-refined mesh, F5, with intersected-element size S = 0.103, is less than F2 with S =

0.069, but the number of computational nodes for F5 is 1.5 times higher than for F2. Another example is between two refined meshes F4 and F7; the refinement scheme for F7 is3 × 3subdividing and for F4 it is4 × 4subdividing. The size of the smallest intersected element for these meshes is the same S = 0.034, but the total number of nodes after refinement for F7 is 1.5 times higher than for F4, while the average difference between the velocity values for these meshes with 1% error is close.

We also sampled the velocity values over the line x = 15.177 mm that crosses the inter-meshing region at a screw orientation ofθ= 112.5○, see Figure 2.16. The trend is similar as in Figure 2.15, except for some oscillations found, and especially clear given the large difference with the BF result in the screw-screw gap region. In the channel the error is about 6-10 % but in the intermeshing gap region the error is large. In this region the drag flow (inydirection) from each screw is in opposite direction and the pressure driven flow dominates the resulting velocity. Apparently this pressure flow is miscalculated in the

−0.80 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 5 10 15 20 25 θ=112.5 x=15.177 mm

y/Rb

V

/V

p

Figure 2.16: The non-dimensional tangential velocity (_VV

p) as a function of

y

Rb over the line

x= 15.177mm for the screw positionθ= 112.5○, using different meshes and

different refinement scheme and comparison with boundary fitted result (BF).

FDM approach.

The total number of nodes can be a good criterion for the size of the memory needed for these computations. With regards to CPU time, in all cases, the refined fictitious domain method is more expensive than the non-refined fictitious domain method. For example for mesh F5 the CPU time with 213257 nodes is 5 times smaller than that for the refined solution F3 with total a number of 201285 nodes. The reason for this is that the refinement procedure takes time.

2.4.4 XFEM results

We solved the same problem with XFEM, using five meshes X1, X2, X3, X4 and X5 summarized in Table 2.4. The convergence of the tangential velocity over the line y =

0,_x < 0 for screw configurationθ = 0○ is shown in Figure 2.17 and also compared to the best result obtained with the finest refined FDM mesh F8. The average error over this line for X5 and BF is less than 2%. Therefore the XFEM predicts the velocity in the channel much better than the fictitious domain method. This is confirmed by comparing the results in the intermeshing region, see Figure 2.18. Also here XFEM shows much more accurate results than the FDM and it predicts the high value of the velocity in ydirection, the same as the BF result. The CPU time for the XFEM is higher than the non-refined fictitious domain method for the same number of elements.

To analyze the large difference between XFEM and FDM results, mainly in the gap region, the pressure values over the intermeshing crossing line are shown for F8, X5 and BF in Figure 2.19. A high-pressure gradient is present in the gap and it causes a local flux. As

Table 2.4: Meshes used for simulations of XFEM.

nelem1 nnodes2 S3 CPU time (s)

X1 1544 6527 0.689 46

X2 5348 22089 0.345 97

X3 13482 55103 0.209 195

X4 29791 120921 0.138 328

X5 52726 213257 0.103 634

1_{Number of elements for reference mesh} 2_{Number of nodes for reference mesh}

3_{Size of intersected elements, S=(circumference of} screw)/(number of intersected elements), circumfer-ence of screw≈ 81.67mm −1 −0.95 −0.9 −0.85 −0.8 −0.75 −0.7 0 0.2 0.4 0.6 0.8 1 y=0 θ=0

x/Rb

V

/V

p

Figure 2.17: The non-dimensional tangential velocity (_VV

p) as a function of

x

Rb over the line

y= 0mm at the left side screw for the screw positionθ= 0○, using different XFEM meshes and comparison with boundary fitted result(BF) and FDM with high-order refinement (F8).

−0.80 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 4 8 12 16 20 23 BF X2 X4 X3 X5 F8 X1 θ=112.5 x=15.177 mm

y/Rb

V

/V

p

Figure 2.18: The non-dimensional tangential velocity (_VV

p) as a function of

y

Rb over the

linex= 15.177mm for the screws positionθ= 112.5○, using different XFEM

meshes and comparison with boundary fitted result (BF) and FDM with high-order refinement (F8).

can be seen in Figure 2.19, this pressure gradient is for the BF and XFEM about 3 times higher than for F8. To verify the effect of the gap region on the accuracy, we increase the gap size between the screws and between screw and barrel by increasing the radius of the barrel to 17 mm and the centerline distance to 28 mm. The results for pressure and velocity in the gap region over the line x = 15.177mm are shown in Figure 2.20. Indeed the differences between XFEM and FDM become much smaller now.

Figures 2.21, 2.22 and 2.23 finally show the contour plots of the velocity in x direction, of the velocity in y direction and the pressure, respectively. Good qualitative agreement found between the XFEM and the BF simulations results. As expected FDM shows larger deviations from the BF results near the high shear-rate regions (gap regions), since the standard finite element integration on intersected elements by the rigid body, cannot pre-dict the jumps in field variables. For the TSE, capturing the jump in the pressure over the gap regions is important to predict the high local value of the velocities and stresses there.

2.5

Conclusions

The objective of this chapter is to compare modified techniques based on mortar ele-ments, the fictitious domain method and the extended finite element method, to simulate the fluid flow inside complex geometries such as those in twin-screw extruders. Using non-conforming mesh refinement, we can simulate the geometries with moving internal

−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 −10 −5 0 5 10 15 F8 BF X5 θ=112.5 x=15.177 mm p y/Rb

Figure 2.19: The relative pressure p (N/mm2) as a function of _Ry

b over the linex= 15.177

mm for the screw positionθ = 112.5○, using refinement mesh F8, X5 and comparison with boundary fitted result (BF).

−0.80 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 1 F X V /V p y/Rb x=15.177 mm θ=112.5 (a) −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 F X y/Rb θ=112.5 x=15.177 mm p (b)

Figure 2.20: Comparison of XFEM (X) and refined FDM (F) meshes for large gap screws for the screw positionθ = 112.5○ over the line x = 15.177mm. (a)The non-dimensional tangential velocity (_VV

p) as a function of

y

Rb. (b)The relative

(a) (b)

(c)

Figure 2.21: Contour color plot of velocity at x direction for screw orientationθ = 0○: (a) BF, (b) X5 and (c) F8.

(a) (b)

(c)

Figure 2.22: Contour color plot of velocity at y direction for screw orientationθ = 157.5○: (a) BF, (b) X5 and (c) F8.

(a) (b)

(c)

Figure 2.23: Contour color plot of pressure for screw orientationθ = 0.0○: (a) BF, (b) X5 and (c) F8.

boundaries just with one structured reference mesh, and refine this mesh locally accord-ing to the position of the movaccord-ing parts. Ensuraccord-ing continuity for field variables in the non-conforming elements is en-forced by using a Lagrangian multiplier. Local mesh refinement increases the accuracy of the velocity field compared to the classical fictitious domain method. Compared to the classical finite element fictitious domain method, local mesh refinement is a much more efficient method. However in more complex geometries, the accuracy of non-conforming mesh refinement is not sufficient anymore. To increase the accuracy of the computed velocity field in e.g. the gap regions of twin-screw extruders, we applied XFEM with some modifications regarding to the geometrical complexity and to the fluid properties. In XFEM the finite element shape functions are extended by using virtual degrees of freedom for the description of the discontinuities around the interface. Then, intersected elements are integrated only on the fluid domain part, and degrees of freedom associated to elements fully inside the rigid body are removed from the sys-tem. The no-slip boundary condition is imposed by using constraints implemented with Lagrangian multipliers. The accuracy and convergence of this scheme has been verified by comparing its results with those of boundary fitted mesh problems. The results are also compared with those based on local mesh refinement in fictitious domain methods. The XFEM method was much more accurate especially in the narrow gap regions, where FDM prove to be inaccurate.