CFD SIMULATION OF A CO-ROTATING TWIN-SCREW EXTRUDER: VALIDATION OF

A RHEOLOGICAL MODEL FOR A STARCH-BASED DOUGH FOR SNACK FOOD

Giorgia Tagliavini(a), Federico Solari(b), Roberto Montanari(c)

(a)
CIPACK Interdepartmental Center, University of Parma – Parco Area delle Scienze, 43124 Parma (Italy)
(b), (c)
Department of Industrial Engineering, University of Parma – Parco Area delle Scienze, 43124 Parma (Italy)

(a)
giorgia.tagliavini@unipr.it, (b)federico.solari@unipr.it, (c)roberto.montanari@unipr.it

ABSTRACT Extrusion is a technique employed for quite some time

The extrusion of starch-based products has been a especially in the food industry. It consists in the
matter of interest, especially for the pasta and snack conveying of high viscosity fluids by increasing the
food production. In recent years, twin-screw extruders pressure from the hopper towards the extruder die.
for snack food were studied from both structural and Extruders can be divided in three key zones: a feeding
fluid dynamics viewpoints. zone, a plasticization zone and a metering zone (as
This project started from the rheological shown in Figure 1).
characterization of a starch-based dough (corn 34wt%,
tapioca 32wt%), comparing viscosity values acquired in
laboratory with different theoretical models found in
literature.
A CFD simulation recreating the simple case of a fluid
flow between two parallel plates was carried out in
order to validate the former comparison. After the Figure 1: Simplified scheme of a co-rotating
rheological characterization was completed, the second intermeshing twin-screw extruder (modified from:
phase of this work covered a 3D CFD simulation of the http://plastictechnologies.blogspot.de).
first part of the twin-screw extruder (feeding zone).
The objective was to find a suitable model for The extruder may also be a single- or twin-screw: in the
describing the dough rheological behavior and the latter, the two parallel screws can co-rotate or counter-
operating conditions of a co-rotating intermeshing twin- rotate, according to the pressure needed (Messori 2005).
screw extruder. A co-rotating twin-screw extruder provides a more
Once the model would be design, it would allow to accurate regulation on the process and a higher range of
investigate several working conditions and different shear stresses (and of shear rate as a result).
screws geometries of the machine, predicting the In the food industry the main goal is to achieve a
evolution of the product rheological properties. smooth and homogeneous products, reason why twin-
screw extruders are preferred over single-screw ones.
Keywords: twin-screw extruders, starch-based dough, Nevertheless, inside twin-screw extruders, it is more
computational fluid-dynamics, rheological model complicated to evaluate the fluid dynamics behavior
and the share rate distributions. For this reason, fluid
1. INTRODUCTION dynamics simulations are essential especially in the
In the food industry, the extrusion of starch-based design phase.
doughs has been widely studied. In the last decade, twin-screw extruders for snack food
One of these doughs main features is the non- were studied from both structural (with finite element
Newtonian behavior, which is characterized by changes method) and fluid-dynamics (with CFD simulations)
in the rheological properties depending on the perspectives (Yamsaengsung et al., 2010; Fabbri et al.,
undergoing process (Pessini, 2001). These rheological 2012; Cubeddu et al., 2014; Bertrand et al., 2003).
changes are affected by temperature, shear stress and Computational fluid dynamics (CFD) has been gaining
shear rate (Lagarrigue & Alvarez, 2001), moisture highlights in recent years for predicting fluid dynamic
content (Rao et al., 1997) and other occurring behavior for a wide range of industrial processes.
substances inside the dough (Xei et al., 2009; Jamilah et The objective of this work is the rheological
al., 2009). characterization of a starch dough for snack food and
All these factors are thus crucial in terms of obtaining the design of a CFD model with ANSYS FLUENT 17.0
the desired product. of a twin-screw extruder first section (feeding zone).

Proceedings of the International Food Operations and Processing Simulation Workshop 32

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.
  u j ui 
2. MATERIALS AND METHODS D    (2)
 xi x j 
 
2.1. Data collection where u and x are the velocity and the coordinate
The first step of this work started with laboratory respectively.
measurements, collecting the viscosity data of the The situation is different when it comes to non-
product at known shear rates. Newtonian fluids. In this case, there is a dependence
The instrumentation used for the data collection between the viscosity  (non-Newtonian viscosity) and
consisted in a Brookfield Engineering RST Controlled
the rate-of-deformation tensor D (Equation 3):
Stress Rheometer (RST-CC Coaxial Cylinder model).
The rheometer and the spindle technical specifications
are listed in Table 1 and Table 2.   ( D )D (3).

Table 1: RST-CC rheometer technical specifications. In relation to this study case, the viscosity  depends
Viscosity Range [Pa∙s] 0.00005 – 5.41M only on the shear rate  , which is related to the tensor
Speed [rpm] 0.01 – 1.3K
Max. Torque [mNm] 100 D according to the following equation:
Torque Res. [µNm] 0.15
1
  D :D (4).
Table 2: CCT-25 spindle technical specifications. 2
Viscosity Range [Pa∙s] 0.002 – 177K
Shear Rate [s-1] 0.013 – 1.67K For simplification purposes, the temperature effects on
Max. Shear Stress [Pa] 2.28K viscosity were omitted at this first stages.
Sample Volume [mL] 16.8 The related terms will, therefore, not be shown in the
following equations.
The viscosity values and the associated shear rates An extensive bibliographic research revealed that the
collected during the laboratory campaign are shown in most commonly used models for describing starch-
Table 3. based fluids are the non-Newtonian Power Law, the
Carreau and the Cross model (Emin & Schuchmann,
Table 3: Collected data from laboratory tests. 2012).
The governing equations of the models mentioned
Shear Rate [s-1] Viscosity [Pa∙s]
above are described below, in accordance with the
4.5 0.171 formulation presented in the ANSYS FLUENT Theory
7.47 0.167 Guide (Ansys, Inc., 2009).
12.51 0.158 2.2.1. Power Law for non-Newtonian Viscosity
20.88 0.137 A non-Newtonian flow is modeled with Equation 5,
40.77 0.115 according to the power law for the non-Newtonian
viscosity:
58.05 0.104
  K n1 (5)
2.2. Rheological models
The next step was the rheological characterization of the
starch (corn 34wt%, tapioca 32wt%) dough. where K is the consistency index, n is the power law
As mentioned before, one of this doughs main features index, a measure of the average viscosity of the fluid
is the non-Newtonian behavior, which implies changes and a measure of the fluid deviation from Newtonian
in the rheological properties, depending on the extrusion behavior respectively, which determines the fluid class:
conditions.
For Newtonian fluids, the shear stress is described by  n  1  Newtonian fluid;
Equation 1, as a function of the rate-of-deformation  n  1  shear-thickening (dilatant fluid);
tensor D :  n  1  shear-thinning (pseudo-plastic).

2.2.2. The Carreau Model for pseudo-plastics

  D (1) The power law model described in Equation 5 results in
a fluid viscosity that varies with shear rate.
where  is the viscosity, which is independent of D . The Carreau model attempts to describe a wide range of
fluids thanks to a curve-fit, which synthetizes functions
D is defined by Equation 2 as: for both Newtonian and shear-thinning non-Newtonian
laws. In the Carreau model, the viscosity is:

Proceedings of the International Food Operations and Processing Simulation Workshop 33

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.
     ( 0    )[ 1   2 2 ] ( n1 ) / 2 (6)

where the parameters n ,  ,  0 and   depend upon

the fluid.  is the time constant, n is the power law
index,  0 and   are respectively the zero- and the
infinite- shear viscosities. The viscosity is limited by
 0 for   0 and by   for    .

2.2.3. The Cross Model Figure 2: Models comparison from the Microsoft Office
The Cross model for viscosity is: Excel elaboration.

0 3. CFD SIMULATION FOR

 (7).
1   
1n RHEOLOGICAL MODELS VALIDATION
A CFD simulation was helpful to identify the most
suitable rheological model among the ones mentioned
 0 is the zero-shear viscosity,  is the natural time previously. A 2D geometry was implemented in
(i.e., inverse of the shear rate at which the fluid changes ANSYS FLUENT, recreating the simple case of a fluid
from Newtonian to power-law behavior) and n is the flow between two fixed, parallel plates.
power law index. Figure 3 shows a detail of the 2D fluid domain used in
the simulation.
2.3. Data processing and analysis
The data were processed with Microsoft Office Excel
2013.
The experimental viscosities were used to derive the
coefficients of the rheological models by minimizing
Figure 3: 2D computational domain with inlet velocity.
the mean squared error between the experimental results
and those obtained using Equations 5, 6 and 7.
The length, in particular, was extended at 5 m to get a
Then the coefficient of determination R2 was estimated.
fully-developed flow inside the duct.
This was a fundamental step for the continuation of the
The domain sizes are presented in Table 6, while Table
study, since the coefficients values are required by the
7 lists the mesh specification.
CFD simulation software for the fluid characterization.
The results are shown in Tables 4 and 5 and Figure 3.
Table 6: Domain sizes.
Table 4: Viscosity values. Size [m]
Height (D) 0.1
Viscosity [Pa∙s] Length (L) 5
Lab tests Power Law Carreau Cross
0.171 0.179 0.172 0.173 The duct was decomposed with ANSYS MESHING,
0.167 0.163 0.166 0.165 according to the Finite Element Method.
0.158 0.148 0.155 0.154 Table 7 lists the mesh specifications.
0.137 0.134 0.139 0.140
0.115 0.118 0.115 0.116 Table 7: Mesh specifications.
0.104 0.110 0.103 0.102 Element Type Quadrilaterals
Num. Elements 5240
Num. Nodes 5523
Table 5: Coefficients and R2 values.
Min. Size 0.0025 m
Power Law Carreau Cross
Max. Size 0.498 m
K 0.239
n Max. Face Size 0.249 m
0.809 0.635 0.129
Wall Edge Sizing Num. of division: 250
 0.083 0.014
Inlet/Outlet Edge Sizing Num. of division: 20
0 0.176 0.189
 0.011 Three different simulations (one for each rheological
R2 0.936 0.997 0.991 model) was carried out in ANSYS FLUENT 17.0 with
stationary operating conditions.
The inlet velocity was fixed at 2.5 m/s and the
atmospheric pressure was set at the outlet.

Proceedings of the International Food Operations and Processing Simulation Workshop 34

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.
For the fluid characterization, the software requires the
coefficients of the chosen model and the density. The
models’ coefficients and a density of 588.94 kg/m3 were
implemented in the simulator.
The fluid flow was set as laminar, due to the high
viscosity of the product.

3.1. Results
The aim of this phase was finding the model that best
fits the experimental data.
First the viscosity values were observed with the
contour using CFD-Post software. Figure 6: Comparison between the viscosity simulation
The results are shown in Figures 4 and 5. and experimental data as a function of shear rate.

The values reveal that the most suitable model was the
Carreau one.
The non-Newtonian Power Law marks the experimental
data quite well, while the Cross model was the less
appropriate.

4. CFD SIMULATION OF THE TWIN-

Figure 4: Dynamic viscosity contour comparison
SCREW EXTRUDER
between rheological models.
The second phase of the study began with 3D modeling
of the twin-screw extruder (feeding zone section).
The geometry reproduced the layout of the extruder’s
first section of the pilot plant, on which the
experimental tests will be carried out to collect a further
amount of data for future validation.
Figure 7 shows the size of the model and Figure 8
depicts the 3D geometry.

Figure 5: Comparative graph with velocity and viscosity

values for different rheological model, depending upon
duct height.

The profiles shown in Figure 5 are dependent upon the

domain height and reveal how well each model
reproduce the fluid behavior in terms of velocity and
viscosity.
Afterwards, the viscosity values were plotted ad a Figure 7: Extruder sizes in (m).
function of the shear rate together with the experimental
data. The shear rate was solved taking the derivative of
the velocity in ∂y.
This approach allows to see which model implemented
in the software would be more accurate.
Finally, the viscosities were put together in the same
graph as a function of ∂v/∂y [s-1].
Figure 6 illustrates the results.

Figure 8: 3D geometry of the twin-screw extruder

implemented in ANSYS FLUENT 17.0.

Proceedings of the International Food Operations and Processing Simulation Workshop 35

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.
The designed geometry defining the computational fluid The atmospheric pressure was set at the inlet, so that the
volume was decomposed with ANSYS MESHING, fluid flow would be influence solely by the screw
according to the Finite Element Method. rotation. A sharp restriction at the outlet cross-section
Special attention was paid to the creation of the mesh, (not shown in Figure 8) was imposed to reproduce the
in order to have accurate information about fluid pressure drop caused by the machine consecutive
dynamic behavior in near-wall regions and in leakage sectors.
areas. Actually, there are narrow gaps between the To impose the angular velocity of 370 rpm a UDF in C
screws and between the screw threads and the barrel. language was created and compiled inside ANSYS
Hence, an inflated mesh is required in those areas to FLUENT.
prevent the lack of details during the calculation. To recreate the screws rotation the dynamic mesh
Furthermore, a dedicated mesh was generated for the option was chosen. This method allows to subdivide the
boundary layer in connection with the solid parts (barrel mesh motion zones, remeshing zones and static zones.
and screws walls). The remeshing zone, precisely, are the areas where the
Table 8 and Figure 9 show the mesh structure and mesh is rebuilt following the solid parts rotation step by
specifications. step during the calculation.
In the present case, before starting to create the dynamic
mesh zones, the fluid domain was split into several
Table 8: Extruder mesh specifications.
parts. Specifically, the boundary layers related to the
Element type Volume: tetrahedrons screws and the barrel walls was separated from the rest
Boundary: quadrilaterals of the fluid domain.
Num. Elements 2721599 In this way the fluid volume was divided into 4 regions.
Num. Nodes 9259878 The areas created by the dynamic mesh are illustrated in
Min. Size 8∙10-5 m Table 9.
Max. Size 1∙10-3 m The term “rigid body” specifies the moving zones,
Max. Face Size 1∙10-3 m “deforming” the zones affected by the remeshing
Inflation (boundary layer) Total thickness process and “stationary” the fixed zones, not subjected
Num. of layers: 4 to remeshing.
Max thickness: 1∙10-4 m
Threads sizing 2∙10-4 m Table 9: Dynamic mesh zones.
Screw sizing 6∙10-4 m Screws Rigid body
Barrel sizing 8∙10-4 m Screws boundary layer Rigid body
Edge sizing 4∙10-4 m Barrel wall Stationary
Barrel wall boundary layer Stationary
Inlet Stationary
Inlet boundary layer Stationary
Outlet Stationary
Outlet boundary layer Stationary
Extruder interior Deforming

Therefore, the dynamic mesh allows the various

boundary layers to move together with its
corresponding wall, so as to facilitate the remeshing
operation and confine it to tetrahedral elements, which
are easier to manage.
The simulation time was chosen in accordance with the
real dough’s residence time inside the machine.

4.1. Results
Once the simulations ended, velocity and viscosity
Figure 9: Mesh structure with details of near-wall values were observed.
region and leakage areas. The velocity vectors (Figure 10) show how the product
is conveyed only by the screws’ rotation. This agrees
The simulation was carried out in a transient state to with the real case and with the boundary conditions
stress the influence of the screw rotation on the flow entered in the software.
behavior. The fluid parameters are those used in the 2D
simulation, i.e. the Carreau model and the laminar flow
regime.

Proceedings of the International Food Operations and Processing Simulation Workshop 36

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.
 5. CONCLUSIONS
This work has both the objective to find a rheological
model for a corn and tapioca dough used for snack food
and design a 3D computational model, which could
recreate the operating conditions of a co-rotating twin-
screw extruder.
A set of viscosity values was obtained experimentally at
specific shear rates.
The experimental data were then compared with the
theoretical viscosity from non-Newtonian Power Law,
Carreau and Cross models, which equations are often
used for describing non-Newtonian fluid behavior.
To find which model was the most suitable, a simple 2D
geometry of two parallel flat planes was created. Three
simulations were carried out with this geometry and it
turned out that the Carreau model was the one that best
reproduced the trend of the experimental values.
Furthermore, the operating condition of the twin-screw
extruder were simulated using the Carreau model for
describing the flow behavior.
Figure 10: Velocity module and vectors. The fluid flows The viscosity values from this simulation were plotted
along the z-positive direction. as a function of chosen sampling sections.
The next step of this project will be the CFD model
validation thanks to an experimental campaign on a
scaled pilot plant, which is now under construction.
This campaign will contemplate the viscosity evaluation
in each sampling section identified during the
computational analysis.
Once the validation will be done, the CFD model could
be used for studying different screws geometries to find
which configuration provides the best results in terms of
product quality.

ACKNOWLEDGMENTS
This work is partly supported by “POR-FESR 2014-
Figure 11: Viscosity contour. 2020, Emilia-Romagna Region, Asse 1: Ricerca e
Innovazione - Bando per progetti di ricerca industriale
Observing values of the dynamic viscosity in Figure 11, strategica rivolti agli ambiti prioritari della Strategia di
it is evident how they agree with the ones found in Specializzazione Intelligente - Progetto: Nuovi
Table 4, for the laboratory results. paradigmi per la progettazione, costruzione ed il
In Figure 11 the sampling sections are shown, from funzionamento di macchine e impianti per l'industria
which the viscosity values were calculated as reported alimentare”.
in the graph in Figure 12.
REFERENCES

Ansys, Inc, (2009), Ansys Fluent Theory Guide.

Chapter 8, Paragraph 4.5.
Bertrand, F., Thibault, F., Delamare, L., & Tanguy, P.
A. (2003). Adaptive finite element simulations of
fluid flow in twin-screw extruders. Computers &
Chemical, 27, 491–500.
Cubeddu, A., Rauh, C., & Delgado, A. (2014). 3D
Thermo-Fluid Dynamic Simulations of High-
Speed-Extruded Starch Based Products. Open
Journal of Fluid Dynamics, 04(01), 103–114.
Figure 12: Viscosity values as function of the sampling Emin, M. a., & Schuchmann, H. P. (2012). Analysis of
sections. the dispersive mixing efficiency in a twin-screw
extrusion processing of starch based matrix.
Journal of Food Engineering, 115(1), 132–143.

Proceedings of the International Food Operations and Processing Simulation Workshop 37

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.
Fabbri, A., Angioloni, A., Di Stefano, A., Fava, E., Roberto MONTANARI is a Full professor of
Guarnieri, A., & Lorenzini, G. (2012). Preliminary Mechanical Plants at the University of Parma. He
Investigation of Pasta Extrusion Process: graduated (with distinction) in 1999 in Mechanical
Rheological Characterization of Semolina Dough. Engineering at the University of Parma. His
Journal of Agricultural Engineering, 38(2), 21–24. research activities mainly concern equipment
Jamilah, B., Mohamed, A., Abbas, K. a., Abdul maintenance, power plants, food plants, logistics,
Rahman, R., Karim, R., & Hashim, D. M. (2009). supply chain management, supply chain modelling
Protein-starch interaction and their effect on and simulation, inventory management. He has
thermal and rheological characteristics of a food published his research in approx. 70 papers, which
system: A review. Journal of Food, Agriculture appear in qualified international journals and
and Environment, 7(2), 169–174. conferences. He acts, as a referee, for several
Lagarrigue, S., & Alvarez, G. (2001). The rheology of scientific journals, is editorial board member of
starch dispersions at high temperatures and high two international scientific journals and editor of a
shear rates: A review. Journal of Food scientific journal.
Engineering, 50(4), 189–202.
Messori, M. (2005). Tecnologie di Trasformazione dei
Materiali Polimerici: Estrusione ed applicazione
dell’estrusione. Stampaggio ad iniezione
Soffiaggio di corpi cavi, 8(I).
Peressini, D. (2001). Evaluation methods of dough
rheological properties. Tecnica Molitoria.
Rao, M. A., Okechukwu, P. E., Da Silva, P. M. S., &
Oliveira, J. C. (1997). Rheological behavior of
heated starch dispersions in excess water: role of
starch granule. Carbohydrate Polymers, 33(4),
273–283.
Xie, F., Yu, L., Su, B., Liu, P., Wang, J., Liu, H., &
Chen, L. (2009). Rheological properties of
starches with different amylose/amylopectin ratios.
Journal of Cereal Science, 49(3), 371–377.
Yamsaengsung, R., Noomuang, C., Law, H., Law, B.,
Law, B., Law, C., Extruder, a S. S. (2010). Finite
Element Modeling for the Design of a Single-
Screw Extruder for Starch-Based Snack Products.
Engineering, III, 8–11.

AUTHORS BIOGRAPHY

Giorgia TAGLIAVINI is a Research Assistant at the

Department of Industrial Engineering at the
University of Parma, where she graduated in
Mechanical Engineering in March 2015. Her
research activities focus mainly on CFD
simulation and process modeling of Industrial and
Food Engineering processes.

Federico SOLARI is a PhD in Industrial Engineering

at the University of Parma, from where he got a
Master degree in Food Industry Mechanical
Engineering, discussing a thesis related to the
design of a plant for the volatile compounds
extraction. He achieved his PhD with a thesis
entitled “Advanced approach to the design of
industrial plants by means of computational fluid
dynamics”. He attended several conferences
related to food process, modelling and simulation.
He published several papers on the same topics on
international journal and conferences.

Proceedings of the International Food Operations and Processing Simulation Workshop 38

978-88-97999-83-6; Bruzzone, Longo, Piera and Vignali Eds.