PDF
Pasta Extrusion
Introduction
The pasta extrusion process consists of several steps. After mixing semolina and water, the hydrated semolina is discharged into an extrusion barrel. In the first part of the extruder, the compression and homogenization sections, the dough is in a granular state. In the next section, the metering zone of the extruder, the material is fully melted and the domain is fully filled with the melt.
This example shows how to simulate the non-isothermal flow of the dough in the metering zone of the pasta extruder taking into account temperature dependent material properties of the hydrated semolina dough.
Model Definition
This model is inspired by Ref. 1, and it consist of a screw that feeds the dough through the extrusion bell toward the extrusion holes, where the spaghetti is made. The screw consist of 10 turns. In the bell, the dough is pressed through a screen with 1-mm holes and further into two outlets with 7 spaghetti strands in each one. The material swelling after the dough exits the extrusion holes is neglected in this model. The geometry is shown in Figure 1.
Since the viscosity of the melt is high, the Reynolds number, which describes the ratio of inertia to viscous forces, is low. At low Reynolds numbers, viscous forces dominate over inertia forces. Thus, the latter may be neglected in the Navier–Stokes equations and the Creeping Flow interface can be used.
The dough exhibits shear thinning (pseudoplastic) behavior. The power law model is a common choice for modeling the flow of a dough. To avoid an unphysical infinite viscosity at zero shear rate, COMSOL Multiphysics implements the Power law model as,
(1)
where m is the fluid consistency coefficient, n is flow behavior index, denotes a reference shear rate, and is a lower shear rate limit.
Figure 1: Geometry of a pasta screw extruder similar to the one presented in Ref. 1.
Because of viscous heating inside the barrel, the temperature variations during the extrusion operation can be significant. The viscosity of the dough decreases with temperature. To adequately describe the flow, temperature dependent material properties should be used. Rheology experiments (Ref. 1) show that the viscosity is also moisture dependent. The viscosity of the dough decreases with increasing moisture content. As moisture content rises, the influence of the temperature on the viscosity decreases. In this model, the moisture content is assumed to be constant. Temperature dependent coefficients for the non-Newtonian power law model are obtained by interpolation from a set of measured viscosities for a range of temperatures and 30% moisture content.
To solve the equations for conservation of energy, mass and momentum in the fluid domain, the Nonisothermal Flow, Laminar Flow interface is used. The interface contains a coupling between the Creeping Flow and the Heat Transfer in Fluids interfaces, and accounts for viscous dissipation.
The screw is encapsulated in a rotating domain with an angular velocity of 20 rpm. The model utilizes the frozen rotor approach to compute a steady-state flow velocity and a temperature field. A frozen rotor analysis is a memory- and time-efficient steady-state approximation. In a sense, this approach can be described as freezing the motion of the moving part in a given position and then observing the resulting flow field with the rotor in that fixed position.
The inlet is set as an Open boundary with a boundary stress of 2 MPa. At the outlet, the pressure is set to 0 Pa. For the heat-transfer equation, a heat flux with a heat-transfer coefficient of 50 W/(m2·K) and an external temperature of 45°C is set on all outer walls of the extruder. The inlet has a fixed temperature of 45°C, and at the outlet boundaries Outflow conditions are applied.
Results and Discussion
The results of the frozen rotor simulation are shown in Figure 2 through Figure 5.
Figure 4 shows the temperature distribution. The temperature increases due to viscous heat dissipation. The viscous heating is generated in the areas with a high shear rate (Figure 5), near the wall of the barrel and inside the extrusion dies. The heat generated near the barrel wall is convected in the axial direction in a helical path along the screw.
The apparent viscosity decreases with increasing temperature (Figure 4). High values of the viscosity in the low-shear-rate zones near the rotating screw lead to poor mixing of the dough.
The velocity plot (Figure 2) shows an uneven velocity distribution in the extrusion channels. This results in different mass flow rates through the extrusion holes and may lead to varying quality of the final product.
Figure 2: Velocity profile in the screw extruder.
Figure 3: Temperature profile in the screw extruder.
Figure 4: Apparent viscosity.
Figure 5: Shear rate.
Reference
1. F. Sarghini, A. Romano, and P. Masi, “Experimental Analysis and Numerical Simulation of Pasta Dough Extrusion Process,” J. Food Eng., vol. 176, pp. 56–70, 2016.
Application Library path: Polymer_Flow_Module/Tutorials/pasta_extrusion
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Fluid Flow>Nonisothermal Flow>Laminar Flow.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select Empty Study.
6
Click  Done.
Study 1
Frozen Rotor
In the Study toolbar, click  Study Steps and choose Stationary>Frozen Rotor.
Geometry 1
Import 1 (imp1)
1
In the Home toolbar, click  Import.
2
In the Settings window for Import, locate the Import section.
3
Click  Browse.
4
Browse to the model’s Application Libraries folder and double-click the file pasta_extrusion.mphbin.
5
Click  Import.
Form Union (fin)
1
In the Model Builder window, under Component 1 (comp1)>Geometry 1 click Form Union (fin).
2
In the Settings window for Form Union/Assembly, locate the Form Union/Assembly section.
3
From the Action list, choose Form an assembly.
4
Select the Create imprints check box.
5
Click  Build Selected.
Component 1 (comp1)
Rotating Domain 1
1
In the Definitions toolbar, click  Moving Mesh and choose Domains>Rotating Domain.
2
Select Domain 26 only.
3
In the Settings window for Rotating Domain, locate the Rotation section.
4
From the Rotation type list, choose Specified rotational velocity.
5
From the Rotational velocity expression list, choose Constant revolutions per time.
6
In the f text field, type 20[RPM].
Next, define the material properties. The physics interface and the chosen fluid model will suggest which material properties are required to solve the model. Therefore, before creating a semolina dough material, choose the Power law model in the Laminar Flow interface.
Laminar Flow (spf)
Fluid Properties 1
1
In the Model Builder window, under Component 1 (comp1)>Laminar Flow (spf) click Fluid Properties 1.
2
In the Settings window for Fluid Properties, locate the Fluid Properties section.
3
Find the Constitutive relation subsection. From the list, choose Inelastic non-Newtonian.
Materials
Material 1 (mat1)
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, locate the Material Contents section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Fluid consistency coefficient
m_pow
m30(T)
Pa·s
Power law
Flow behavior index
n_pow
n30(T)
1
Power law
Heat capacity at constant pressure
Cp
2000
J/(kg·K)
Basic
Density
rho
1200
kg/m³
Basic
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
0.242
W/(m·K)
Basic
Define interpolation functions for the temperature-dependent power-law coefficients.
Interpolation 1 (int1)
1
In the Model Builder window, expand the Material 1 (mat1) node.
2
Right-click Component 1 (comp1)>Materials>Material 1 (mat1)>Power law (PowerLaw) and choose Functions>Interpolation.
3
In the Settings window for Interpolation, locate the Definition section.
4
In the Function name text field, type m30.
5
In the table, enter the following settings:
t
f(t)
35
109341
40
119257
45
81950
50
41705
55
33043
60
34739
80
33643
6
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
m30
Pa*s
7
In the Argument table, enter the following settings:
Argument
Unit
t
degC
Interpolation 2 (m2)
1
Right-click Interpolation 1 (m30) and choose Duplicate.
2
In the Settings window for Interpolation, locate the Definition section.
3
In the Function name text field, type n30.
4
In the table, enter the following settings:
t
f(t)
35
0.3984
40
0.3365
45
0.3586
50
0.4510
55
0.4540
60
0.4501
80
0.2910
5
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
n30
1
Pasta 30% Hydration
1
In the Model Builder window, under Component 1 (comp1)>Materials>Material 1 (mat1) click Power law (PowerLaw).
2
In the Settings window for Power Law, locate the Output Properties section.
3
In the table, enter the following settings:
Property
Variable
Expression
Unit
Size
Fluid consistency coefficient
m_pow
m30(T)
Pa·s
1x1
Flow behavior index
n_pow
n30(T)
1
1x1
4
In the Model Builder window, click Material 1 (mat1).
5
In the Settings window for Material, type Pasta 30% Hydration in the Label text field.
6
In the Model Builder window, collapse the Component 1 (comp1)>Materials>Pasta 30% Hydration (mat1)>Power law (PowerLaw) node.
The flow velocity in this model is very slow. Thus, it is safe to neglect the inertial terms of the Navier–Stokes equations and use the Stokes equation for creeping flow.
Laminar Flow (spf)
1
In the Model Builder window, under Component 1 (comp1) click Laminar Flow (spf).
2
In the Settings window for Laminar Flow, locate the Physical Model section.
3
From the Compressibility list, choose Incompressible flow.
4
Select the Neglect inertial term (Stokes flow) check box.
To eliminate unphysical behavior of the Power law model at zero shear rate, the COMSOL Multiphysics implementation introduces a cutoff shear rate value. Modify the default minimum shear rate value according to Ref. 1
Fluid Properties 1
1
In the Model Builder window, under Component 1 (comp1)>Creeping Flow (spf) click Fluid Properties 1.
2
In the Lower shear rate limit text field, type 1[1/s].
Open Boundary 1
1
In the Physics toolbar, click  Boundaries and choose Open Boundary.
2
Select Boundary 11 only.
3
In the Settings window for Open Boundary, locate the Boundary Condition section.
4
In the f0 text field, type 2[MPa].
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Select Boundaries 1658–1671 only.
Heat Transfer in Fluids (ht)
In the Model Builder window, under Component 1 (comp1) click Heat Transfer in Fluids (ht).
Outflow 1
1
In the Physics toolbar, click  Boundaries and choose Outflow.
2
Click the  Select Box button in the Graphics toolbar.
3
Select Boundaries 1658–1671 only.
Temperature 1
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
Select Boundary 11 only.
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type 45[degC].
Heat Flux 1
1
In the Physics toolbar, click  Boundaries and choose Heat Flux.
2
Select Boundaries 1, 2, 5, 8, 58, 60, 62–64, 87, 88, 90, 91, 115, 117–120, 122, 123, 1571–1575, 1577, 1578, and 1580–1585 only.
3
In the Settings window for Heat Flux, locate the Heat Flux section.
4
From the Flux type list, choose Convective heat flux.
5
In the h text field, type 50.
6
In the Text text field, type 45[degC].
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Fine.
4
Click  Build All.
Study 1
In the Home toolbar, click  Compute.
Results
Slice
1
In the Model Builder window, expand the Velocity (spf) node, then click Slice.
2
In the Settings window for Slice, locate the Plane Data section.
3
From the Plane list, choose zx-planes.
4
In the Planes text field, type 1.
5
In the Velocity (spf) toolbar, click  Plot.
To improve the visualization, it is possible to include the metal surfaces of the screw and render it as a scratched steel surface. To facilitate this, it is best to define an explicit selection of the desired surfaces.
Definitions
Metal Surfaces
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, locate the Input Entities section.
3
From the Geometric entity level list, choose Boundary.
4
Click the  Wireframe Rendering button in the Graphics toolbar.
5
Select Boundaries 8, 15, 16, 19, 20, 23, 24, 27, 28, 31, 32, 35, 36, 39, 40, 43, 44, 47, 48, 51, 52, 57, 58, 62, 67, 69, 71, 73, 75, 77, 79, 81, 87–89, 94, 96, 98, 100, 102, 104, 106, 108, 110, 112, 119, 120, 123, 124, 126–128, 130–132, 134–136, 138–140, 142–149, 151–153, 155–157, 159–161, 163–165, 167–169, 171–173, 175–177, 179–181, 183–185, 187–199, 201–203, 205–207, 209–211, 213–215, 217–219, 221–223, 225–227, 229–231, 233–235, 237–239, 241–243, 245–259, 261–263, 265–267, 269–271, 273–275, 277–279, 281–283, 285–287, 289–291, 293–295, 297–299, 301–303, 305–307, 309–311, 313–315, 317–334, 336–338, 340–342, 344–346, 348–350, 352–354, 356–358, 360–362, 364–366, 368–370, 372–374, 376–378, 380–382, 384–386, 388–390, 392–394, 396–414, 416–418, 420–422, 424–426, 428–430, 432–434, 436–438, 440–442, 444–446, 448–450, 452–454, 456–458, 460–462, 464–466, 468–470, 472–474, 476–478, 480–499, 501–503, 505–507, 509–511, 513–515, 517–519, 521–523, 525–527, 529–531, 533–535, 537–539, 541–543, 545–547, 549–551, 553–555, 557–559, 561–563, 565–567, 569–589, 591–593, 595–597, 599–601, 603–605, 607–609, 611–613, 615–617, 619–621, 623–625, 627–629, 631–633, 635–637, 639–641, 643–645, 647–649, 651–653, 655–657, 659–661, 663–684, 686–688, 690–692, 694–696, 698–700, 702–704, 706–708, 710–712, 714–716, 718–720, 722–724, 726–728, 730–732, 734–736, 738–740, 742–744, 746–748, 750–752, 754–756, 758–779, 781–783, 785–787, 789–791, 793–795, 797–799, 801–803, 805–807, 809–811, 813–815, 817–819, 821–823, 825–827, 829–831, 833–835, 837–839, 841–843, 845–847, 849–851, 853–874, 876–878, 880–882, 884–886, 888–890, 892–894, 896–898, 900–902, 904–906, 908–910, 912–914, 916–918, 920–922, 924–926, 928–930, 932–934, 936–938, 940–942, 944–946, 948–969, 971–973, 975–977, 979–981, 983–985, 987–989, 991–993, 995–997, 999–1001, 1003–1005, 1007–1009, 1011–1013, 1015–1017, 1019–1021, 1023–1025, 1027–1029, 1031–1033, 1035–1037, 1039–1041, 1043–1064, 1066–1068, 1070–1072, 1074–1076, 1078–1080, 1082–1084, 1086–1088, 1090–1092, 1094–1096, 1098–1100, 1102–1104, 1106–1108, 1110–1112, 1114–1116, 1118–1120, 1122–1124, 1126–1128, 1130–1132, 1134–1136, 1138–1159, 1161–1163, 1165–1167, 1169–1171, 1173–1175, 1177–1179, 1181–1183, 1185–1187, 1189–1191, 1193–1195, 1197–1199, 1201–1203, 1205–1207, 1209–1211, 1213–1215, 1217–1219, 1221–1223, 1225–1227, 1229–1249, 1251–1253, 1255–1257, 1259–1261, 1263–1265, 1267–1269, 1271–1273, 1275–1277, 1279–1281, 1283–1285, 1287–1289, 1291–1293, 1295–1297, 1299–1301, 1303–1305, 1307–1309, 1311–1313, 1315–1334, 1336–1338, 1340–1342, 1344–1346, 1348–1350, 1352–1354, 1356–1358, 1360–1362, 1364–1366, 1368–1370, 1372–1374, 1376–1378, 1380–1382, 1384–1386, 1388–1390, 1392–1409, 1411–1413, 1415–1417, 1419–1421, 1423–1425, 1427–1429, 1431–1433, 1435–1437, 1439–1441, 1443–1445, 1447–1449, 1451–1453, 1455–1457, 1459–1461, 1463–1479, 1481–1483, 1485–1487, 1489–1491, 1493–1495, 1497–1499, 1501–1503, 1505–1507, 1509–1511, 1513–1515, 1517–1529, 1531–1533, 1535–1537, 1539–1541, 1543–1545, 1547–1549, 1551–1553, 1555–1557, 1559–1569, 1573, 1574, 1582–1586, 1593, 1675, 1678, 1698–1720, and 1731–1752 only.
This is most easily done by copying the selection from the text above and pasting it in the selection area.
6
In the Label text field, type Metal Surfaces.
Next, add a surface plot with the previous selection. Then you can add a material appearance to render it as a scratched steel surface.
Results
Surface 1
1
In the Model Builder window, right-click Velocity (spf) and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type 1.
4
Click to expand the Title section. From the Title type list, choose None.
Selection 1
1
Right-click Surface 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Metal Surfaces.
Material Appearance 1
1
In the Model Builder window, right-click Surface 1 and choose Material Appearance.
2
In the Settings window for Material Appearance, locate the Appearance section.
3
From the Appearance list, choose Custom.
4
From the Material type list, choose Steel (scratched).
Copy the surface plot with the metal rendered screw for use in other plot groups.
Surface 1
Right-click Surface 1 and choose Copy.
Modify the temperature plot so that it looks similar to Figure 3 by deleting the original surface plot and adding a temperature slice instead.
Surface
1
In the Model Builder window, expand the Results>Temperature (ht) node.
2
Right-click Results>Temperature (ht)>Surface and choose Delete.
Slice 1
1
In the Model Builder window, right-click Temperature (ht) and choose Slice.
2
In the Settings window for Slice, locate the Expression section.
3
In the Expression text field, type T.
4
From the Unit list, choose degC.
5
Locate the Plane Data section. From the Plane list, choose zx-planes.
6
In the Planes text field, type 1.
7
Locate the Coloring and Style section. Click  Change Color Table.
8
In the Color Table dialog box, select Thermal>ThermalDark in the tree.
9
Click OK.
Surface 1
Right-click Temperature (ht) and choose Paste Surface.
Add a new plot group and create a slice plot for the visualization of the viscosity, similar to Figure 4. Include the previous surface plot of the steel screw parts.
3D Plot Group 5
In the Home toolbar, click  Add Plot Group and choose 3D Plot Group.
Surface 1
Right-click 3D Plot Group 5 and choose Paste Surface.
3D Plot Group 5
In the Model Builder window, right-click Surface 1 and choose Slice.
Slice 1
1
In the Settings window for Slice, locate the Plane Data section.
2
From the Plane list, choose zx-planes.
3
Locate the Expression section. In the Expression text field, type spf.mu.
4
Locate the Plane Data section. In the Planes text field, type 1.
5
In the 3D Plot Group 5 toolbar, click  Plot.
Viscosity
1
In the Model Builder window, under Results click 3D Plot Group 5.
2
In the Settings window for 3D Plot Group, type Viscosity in the Label text field.
Duplicate the plot group and add a plot for the shear rate.
Shear Rate
1
Right-click Viscosity and choose Duplicate.
2
In the Settings window for 3D Plot Group, type Shear Rate in the Label text field.
Slice 1
1
In the Model Builder window, expand the Shear Rate node, then click Slice 1.
2
In the Settings window for Slice, locate the Expression section.
3
In the Expression text field, type spf.sr.
4
In the Shear Rate toolbar, click  Plot.
5
In the Expression text field, type log10(spf.sr).
6
In the Shear Rate toolbar, click  Plot.