DOI: 10.3303/CET25118004 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Paper Received: 28 May 2025; Revised: 23 June 2025; Accepted: 30 August 2025 
Please cite this article as: Di Lorenzo G., Albanese D., Marra F., 2025, Preliminary Step of the Digital Twin Model Design for High Moisture 
Extrusion of Proteins from Alternative Sources, Chemical Engineering Transactions, 118, 19-24  DOI:10.3303/CET25118004 
  

 CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 118, 2025 

A publication of 

 

The Italian Association 

of Chemical Engineering 

Online at www.cetjournal.it 

Guest Editors: Laura Piazza, Francesco Donsì, Giorgia Spigno, Paolo Ciambelli, Luca di Palma 

Copyright © 2025, AIDIC Servizi S.r.l. 

ISBN 979-12-81206-19-9; ISSN 2283-9216 

Preliminary Step of the Digital Twin Model Design for High 

Moisture Extrusion of Proteins from Alternative Sources 

Gianfranco Di Lorenzo*, Donatella Albanese, Francesco Marra 

Industrial Engineering Department, University of Salerno, Fisciano (SA), Campania 84084, Italy.  

gdilorenzo@unisa.it  

The objective of this study is to develop a deeper understanding of the high-moisture extrusion (HME) process 

for producing meat-like products from alternative proteins and to establish a foundation for the implementation 

of a Digital Twin (DT) framework. Since twin screw extruder (TSE) is largely a black-box system with limited 

temperature sampling and no direct velocity measurements, a mechanistic modelling framework was developed 

as an initial step toward DT implementation. To achieve this, numerical simulations were conducted to analyse 

key transport phenomena and axial velocity distribution, pumping efficiency, and dispersive mixing behaviour, 

within a co-rotating twin screw extruder. A steady-state approximation with frozen rotor analysis was employed 

to model TSE behaviour. The numerical simulations provided insights into the axial velocity distribution along 

the screws, the effect of screw speed on pumping efficiency, and the role of dispersive mixing in the process. 

These findings contribute to a better understanding of the extrusion process and lay the groundwork for future 

advancements in DT applications, aiming at process optimization, improved material structuring, and enhanced 

production efficiency for alternative protein-based meat analogues. 

1. Introduction 

High-moisture extrusion processing has been used for years to transform plant-based proteins into fibrous, 

meat-like structures at 40-80% moisture, with an extended cooling die aiding fiber formation (Hong et al., 2022). 

To support this, the cooling zone must stay below 75 °C (Xiaonan et al., 2024). Moisture, screw, temperature, 

and pressure affect the material, but detecting changes is challenging. Numerical simulation helps reveal protein 

flow behavior (Dashu et al., 2024). With industrial digitalization and the rise of the Internet of Things, the Digital 

Twin concept has emerged as a tool for process management, though its first agrifood application is yet to be 

demonstrated (Verboven et al., 2020). A Digital Twin, simulating transport phenomena, enables live interaction 

between the physical and digital systems, allowing the machine to learn from real-world data to optimize process 

conditions and develop products with desired characteristics (Datta et al., 2022). This study aims to deepen 

understanding of high-moisture extrusion for meat-like products and lay the groundwork for a Digital Twin 

framework by developing a mechanistic modeling approach to simulate transport phenomena and flow behavior. 

2. Material and methods 

2.1 Design of co-rotating twin screw extruder 

A lab-scale HME system was considered for mechanistic modeling framework design. HME system consisted 

of a parallel co-rotating TSE (L/D: 30 with a cross-sectional diameter of 13.4 mm) fitted with a liquid-cooled slot 

die (H x W x L : 44 x 77 x 270 mm) that leads the extrudate into a slab of rectangular section (H x W: 5 x 18 

mm). Along the extruder, nine alternate heating cartridges (considering the power range of 0- 200 W for each 

cartridge), coupled with temperature sensors, allowed the control of the temperature profile along the extruder 

length. The pressure value was monitored at inlet of liquid-cooled slot die.  

Screws of extrusion 

As shown in Figure 1, the twin screw is comprised of the following elements: 

19

mailto:gdilorenzo@unisa.it


• feeding zone composed of n.2 screw elements with 8 mm lead (8/13.4/0.53/R) and n.2 screw elements 

with 20 mm lead and 0.53 mm crest (20/13.4/0.53/R); 

• transport zone composed of n.15 screw elements with 12 mm lead and 0.53 mm crest (12/13.4/0.53/R); 

• mixing zone composed of kneading block with n.3 elements at 0° orientation (0°/13.4/K) and n.5 

elements with 90° orientation (90°/13.4/K); 

• extrusion zone composed of n.2 elements with 12 mm lead and 2 mm crest (12/13.4/2/R). 

 
Figure 1 Screw configuration 

2.2 Mathematical model development 

The description of the problem under consideration must take into account the balances of momentum, energy 

and matter, with the additional term of viscous generation. The complexity of the problem lies in dealing with an 

intrinsically variable-domain system that was described with frozen rotor analysis. On Table 1, Eqs(1, 2 and 3) 

were used to determine the temperature distribution within the material processed along the extruder. Navier-

Stokes and continuity equations (Eqs(4 and 5)) were coupled with Eqs(1, 2 and 3). In addition, gravity could be 

neglected as the Jeffrey-number (Kousemaker et al., 2024) – the ratio of gravity to viscous forces – is low. 

Table 1 Transport phenomena equation used in mathematical model 

Fourier equation ρ Cp u ∙ ∇T+ ∇∙(-k ∇T)=Qp+Qvd (1) 
Navier Stokes 

equation 
ρ (u∙∇)u=-∇P+ μ (∇u+(∇u)T) (4) 

 

Work done by 

pressure changes 
Qp=-

1

ρ

∂ρ

∂T
 T (u ∙ ∇P) (2) 

Continuity 

equation 
ρ ∇u=0 (5) 

 

Viscous 

dissipation effect 
Qvd=τ : ∇u (3)    

 

Boundary conditions 

To determine the flow field generated by the twin screw motion, a rotation velocity was equally applied to both 

screws. Velocity was changed at 100 RPM, 250 RPM and 500 RPM while the extrusion chamber was kept 

relatively fixed. For the flow field, an open boundary condition was applied at the inlet of the extruder chamber, 

while atmospheric pressure was applied at the outlet of the cooled slot die. Furthermore no-slip condition was 

used along the walls of the extruder chamber. For heat transfer equation, inlet boundary condition was applied 

with a temperature equal to 25°C, while assigning temperatures of 25°C, 25°C, 40°C, 60°C, 80°C, 100°C, 120°C, 

130°C, 110°C respectively at each zone from 1 to 9. Temperature at cooling slot die was set at 70°C to have 

product temperature below 75°C (Xiaonan et al., 2024). 

Material properties 

The geometry block was assumed to be the extruder chamber and cooling slot die completely filled with soy 

dough, having 65 % w/w moisture content. This moisture content was previously used (Lin et al., 2000) to 

extrude in a lab-scale HME process. Physical properties like specific heat Cp (J kg-1 K-1), thermal conductivity k 

(W m-1 K-1) and thermal diffusivity α (m2/s) and density ρ (kg/m3) were obtained by literature (Högg and Cornelia, 

2023) as a function of temperature, except for the density. Viscosity, a key property, strongly depends on the 

fluid dynamics and temperature and this was accounted using parameters, power n and consistency coefficient 

m (Pa s), according to ‘Power Law’ model (Dashu et al., 2024). 

Discretization, mesh and convergence 

A second order (P2) velocity and first order (P1) pressure configuration were selected for the discretization of 

the fluid variables to meet the ‘inf-sup’ (LBB) condition (Arnold et al., 1984) and achieve a smooth and high 

accuracy velocity field. Mesh convergence studies were performed to ensure that the results were independent 

of mesh resolution. With this aim, the influence of meshing a fluid flow field and heat transfer solution, coupled 

with rotor moving mesh, was analyzed using the highest screw speed (500 RPM). Maximum element size and 

minimum element size were set by changing the size mesh parameter (SM [m]) using the parametric sweep 

− 

20



analysis. The outward mass flow rate (Eq(6)) and the viscous dissipation (Eq(3)) were calculated and used to 

evaluate the mesh influence. 

ṁ=ρ ∫u dS  (6) 

According to literature (Ali G., 2013), volumetric flow rates fluctuation between 1 to 8 % should be considered 

to be indicating that the mesh dependency is negligible. The results found that the calculated mass flow rate 

changed between 10.596 to 12.971 kg/h. The presented results always corresponded to the highest mesh 

density value shown in Table 2. 

Table 2 Mesh dependency analysis at different number of cells 

Case Numbers of 

degrees of 

freedom [unit] 

Mass flow rate 

[kg/h] 

Error MFR [%] Total viscous 

dissipation [W] 

Error TVD [%] 

Case 1 1.88·106 12.971 / 473 / 

Case 2 1.23·106 12.191 6 478 1 

Case 3 0.88·106 11.164 8.4 482 0.8 

Pumping efficiency 

The ratio between the total throughput and theoretical throughput based on its C-shaped chamber volume is 

called pumping efficiency (Ali G., 2013). This factor multiplied the theoretical throughput gives the total 

throughput. The theoretical throughput Qth was calculated as reported in Eq(7) (Ali G., 2013). 

Qth=2 N Vc ρ  (7) 

The C-shaped chamber volume can be calculated by subtracting the volume of a given length of screw from the 

same length of the empty barrel. The volume V1 of one side of the barrel bore over one flight lead was calculated 

as in Eq(8). 

V1= {(π-
α

2
) R

2
+ (R-

H

2
) √(RH-H

4)} S  (8) 

These geometrical parameters are shown in  Figure 2. 

 

Figure 2 Geometrical parameters for simple model of theoretical output 

α is defined as the overlapping angle (rad) and was calculated as in Eq(9). 

α=2 tan
-1 (

√R
2
-
H

2

4

R-
H

2

)  (9) 

The volume V2 of the screw root over one pinch length was calculated as in Eq(10). 

V2=S π(R-H)2  (10) 

The volume V3, of one screw flight, was calculated as in Eq(11). 

V3=2π [(RH-
H

2

2
) Wf+ (RH

2
-

2

3
H

3) tan(θ)]  (11) 

21



The total volume of C-shaped chamber was calculated as in Eq(12). 

Vc=
V1-V2-iV3

i
  (12) 

Values of the up mentioned variables were calculated and reported in Table 3. 

Table 3 Geometrical parameters and mass flow rate results at different screw speed 

ρ Soy dough density [kg/m3] 1120.00 

i Number of screw flight [-] 2 

R Screw radius [mm] 6.70 

H Channel depth [mm] 1.90 

S Screw lead length [mm] 12.00 

Wf Flight width [mm] 0.53 

θ Flight angle [°] 16.00 

Qth
 Theoretical mass flow rate 

(no leakage) 

100/250/500 

RPM 

[kg/h] 3.13/7.82/15.62 

Dispersive mixing behavior 

Ali G. (2013) mentioned that evaluating the dispersive mixing parameter of a twin screw extruder was a valuable 

measure. The dispersive mixing parameter λ is defined by Eq(8) (Manas-Zloczower et al., 1989). 

λ=
|γ̇|

|γ̇|+|ω|
  (8) 

λ quantifies the elongational and rotational flow which are often thought to be more strongly relevant to 

dispersive mixing compared to shear. It has a value between 0 and 1, where λ = 0 indicates pure rotation while 

λ = 1 indicates pure elongation and λ = 0.5 corresponds to simple shear flow. 

2.3 Numerical solution  

The CFD and heat transfer simulation were carried out using the commercial solver Comsol Multiphysics 6.2. 

Computations were carried out on a PC with Intel® Core™ i7-7700HQ CPU @ 2.80 GHz 2.81 GHz, RAM 16 

GB, NVIDIA GeForce GTX 1050 12228 MB 640 CUDA Core. 

3. Results 

3.1 Axial velocities - effect of screw lead length and tip width 

The screw profile is repeated for each lead length in feeding, transport and extrusion zones. Consequently, only 

axial velocities distribution from 0° to 135 ° planes are plotted. Axial velocities distribution in planes are shown 

at screw speed equal to 500 RPM for feeding zone in Figure 3, for transport zone in Figure 4, and for extrusion 

zone in Figure 5.  Figure 3 for feeding zone shows highest values of axial velocity occurred at 135° plane and 

before the converging section of the intermeshing region. While most negative values are located at 45° plane 

inside the intermeshing region. These results are in according with the literature (Ali G., 2013). The symmetry 

along screws’ inter-axis generates negative values of axial velocities. Similar behaviors are shown in Figure 4 

and Figure 5 suggesting that there is no dependence on the lead length and tip width. However, along transport 

and extrusions zone negative values of velocities are less pronounced, suggesting a favorable lower material 

leakage or backflow. This result potentially affects material residence time and mixing efficiency during extrusion 

process. 

0° (y=-0.02400 m)  45° (y=-0.02645 m) 90° (y=-0.0289 m) 135° (y=-0.03135 m) 

    

Figure 3 Axial velocities distribution at 500 RPM in feeding zone at planes from 0° to 135° (20 mm lead) 

22



 0° (y=-0.1100 m)  45° (y=-0.1115 m) 90° (y=-0.1130 m) 135° (y=-0.1145 m) 

    

Figure 4 Axial velocities distribution 500 RPM in transport zone at planes from 0 to 135° (12 mm lead) 

 0° (y=-0.3655 m) 45° (y=-0.3670 m) 90° (y=-0.3686 m) 135° (y=-0.3701 m) 

    

Figure 5 Axial velocity distribution at 500 RPM in extrusion zone from 0 to 135° planes (12 mm lead) 

3.2 Pumping behaviour – effect of screw speed 

The pumping efficiency was calculated for each screw speed. The results are shown in Table 4. An increase in 

screw speed modified slightly the pumping efficiency (0.77-0.83) suggesting that screw speed affects the 

theoretical flow rate and total flow rate at the same degree according to the literature (Ali G., 2013). This 

suggests that leakage flows and other secondary effects had a negligible impact on throughput. 

Table 4 Calculation of total mass flow rate and theoretical mass flow rate 

Qac Total mass flow rate 100/250/500 RPM [kg/h] 2.40/6.29/12.97 

η Pumping efficiency 100/250/500 RPM / 0.77/0.80/0.83 

3.3 Dispersive mixing behaviour - effect of screw speed 

Quantification of mixing is challenging in extrusion processing. To evaluate the dispersive mixing efficiency the 

dispersive mixing parameter λ was used in this analysis in the mixing zone, its definition is reported in Eq(8). 

|γ|, the magnitude of the rate of strain tensor, and |ω|, the magnitude of the rate of vorticity tensor, were 

calculated by Comsol Multiphysics and results are shown in Table 5. Results show that no influence of screw 

speed on dispersive mixing parameter occurred (0.451-0.454). On the basis of previous mixing parameter 

classification, obtained values suggest a simple shear flow behavior in the mixing zone. 

 Table 5 Dispersive mixing parameter evaluated in the mixing zone at different screw speed 

|γ| 
Magnitude of the rate 

of strain tensor 

100/250/500 RPM [1/s] 60/151/300 

|ω| 
Magnitude of the rate 

of vorticity tensor 

100/250/500 RPM [1/s] 72/180/365 

λ Dispersive mixing 

parameter 

100/250/500 RPM / 0.454/0.456/0.451 

4. Conclusions 

Axial velocity distribution showed that the highest velocities occurred in the intermeshing regions of the screws, 

promoting material transport, while localized negative velocities were observed in specific planes. These 

variations suggest that backflow and recirculation zones are influenced by screw geometry, potentially affecting 

23



material residence time and mixing efficiency with a potential impact on fiber alignment and protein structuring 

during texturization process. Pumping efficiency remained relatively stable (0.77–0.83) across different screw 

speeds, indicating that the total mass flow rate increased proportionally with the theoretical value predicted by 

an analytical model. This suggests that leakage flows and other secondary effects had a negligible impact on 

throughput, confirming the reliability of theoretical predictions for process scaling and optimization. The 

dispersive mixing parameter exhibited minimal variation (0.451–0.456) across different screw speeds, indicating 

that the mixing zone predominantly induced simple shear flow. This behaviour suggests that screw speed does 

not significantly alter dispersive mixing efficiency under the tested conditions. This flow type is beneficial for 

partial protein denaturation and molecular realignment, facilitating the formation of structured protein networks. 

Maintaining a controlled shear environment prevents excessive fragmentation, ensuring the development of 

elongated protein fibers essential for meat analog texture. The study offers useful insights for optimizing 

extrusion parameters, enhancing process control, and improving the consistency of protein structuring in meat 

analog production.  

Nomenclature

ρ – material density, kg/m3Cp – specific heat 

capacity of soy dough, J/(kg·K) 

k – thermal conductivity of soy dough, W/(m·K) 

T – temperature, K 

u – velocity field, m/s 

Qp – work done by pressure changes, W/m3 

Qvd – viscous dissipation in the fluid, W/m3 

P – pressure, Pa 

τ – viscous stress tensor, Pa 

Qth – theoretical mass flow rate, kg/h 

Qac – total mass flow rate, kg/h  

|γ| - magnitude of the rate of strain tensor, 1/s 

|ω| - magnitude of the rate of vorticity tensor, 1/s 

μ – non-newtonian viscosity, Pa·s 

N – screw rotation speed, 1/s 

Vc – C-shaped volume chamber, m3 

η – pumping efficiency, - 

S – screw lead length, m 

R - screw radius, m 

H – channel depth, m 

Wf – flight width, m 

α – overlapping angle, ° 

θ – flight angle,  ° 

y – axis direction, m

References 

Ali G. (2013). Modelling of counter rotating twin screw extrusion, Master Thesis, Mc Master University, 1-145. 

Copyright by Ali Goger, August 2013, Hamilton, Ontario. 

Arnold, D.N., Babuˇska, I., Osborn, J., 1984. Finite element methods: Principles for their selection. Comput. 

Methods Appl. Mech. Eng. 45, 57–96.  

Dashu L., Tianqui Z., Fei G., Chengyi Z., Dongyu S., Yutong G and Min W.(2024). The study of numerical 

simulation and texture of soybean protein based on high moisture extrusion with different screw elements. 

Innovative Food Science and Emerging Technologies, 92, 103560. 

Datta A., Nicolai B., Vitrac O., Verboven P., Erdogdu F., Marra F., Sarghini F., Koh C. (2022). Computer-aided 

food engineering. Nature Food, Nov;3(11):894-904. 

Högg E. and Cornelia R. (2023). Towards a better understanding of Texturization during High-Moisture 

Extrusion. Foods, 12, 2283.  

Hong W., Frans W.J.V., Wei Z., Tomasz P.C., Longteng Z., Birthe M.J., R.L. (2022). Differences in 

physicochemical properties of high-moisture extrudates prepared from soy and pea protein isolates. Food 

Hydrocolloids, 128, 107540. 

Kousemaker, T. M., Druetta, P., Picchioni, F., & Vakis, A. I. (2024). Modelling supercritical CO2 flow in a co-

rotating twin screw extruder using the level-set method. Chemical Engineering Research and Design, 205, 

569-577. 

Lin. S., Huff H.E. and Hsieh F. (2000). Texture and Chemical Characteristics of Soy Protein Meat Analog 

Extruded at High Moisture. Journal of Food Science, Vol.65, No.2.  

Manas-Zloczower I. and Feke D.L.(1989). Int. Polym. Proc., IV, 3-8. 

Verboven, P., Defraeye, T., Datta, A. K. & Nicolai, B. (2020). Digital twins of food process operations: the next 

step for food process models? Curr. Opin. Food Sci. 35, 79–87. 

Xiaonan S., Tianyi Z., Xin Z., Lianzhou J. (2024). High-Moisture Extrusion of Plant Proteins: Fundamentals of 

Texturization and Applications. Annu. Rev. Food Sci. Technol., 15:125-49. 

24


	10dilorenzo.pdf
	Preliminary Step of the Digital Twin Model Design for High Moisture Extrusion of Proteins from Alternative Sources