Extrusion Simulation for the Design of Cereal and Legume Foods

A 1D global twin-screw extrusion model, implemented in numerical software, Ludovic®, was applied to predict extrusion variables and, therefore, to design various starchy products with targeted structure and properties. An experimental database was built with seven starchy food formulations for manufacturing dense and expanded foods made from starches, starch blends, breakfast cereals, pulse crop ingredients such as pea flour, fava bean flour, and fava bean starch concentrated, and wheat flour enriched with wheat bran. This database includes the thermal and physical properties of the formulations at solid and molten states, melt viscosity model, extruder configurations and operating parameters, and extruded foods properties. Using extrusion and viscosity models, melt temperature (T) and specific mechanical energy (SME) were satisfactorily predicted. A sensitivity analysis of variables at die exit was performed on formulation, extruder configuration, and operating parameters, generating the extruder operating charts. Results allowed the establishment of relationships between predicted variables (T, SME, melt viscosity) and product features such as starch and protein structural change, density and cellular structure, and functional properties. The extrusion operating conditions leading to targeted food properties can be assessed from these relationships and also the relationship between extrusion operating parameters and variables provided by simulation.


Introduction
Many starch-based foods, such as texturized ingredients, breakfast cereals, snacks, crackers, pasta, noodles, pet food, and many others, are produced by extrusion processes [1]. Starches and their blends are extruded to gain desirable end-usage properties. Ready-to-eat breakfast cereals are still the market leader in extruded foods. Society is leaning increasingly towards more sustainable and healthy diets by transitioning from animal to plant protein sources and enriching formulations with legume proteins and dietary fibers. Despite their nutritional advantages (high protein content, low glycaemic index) compared to cerealbased extruded foods, legume-based extruded foods are still rare. The growing demand for extruded snacks (4.4% per year in the market) makes the development of healthy extruded snacks even more desirable [2].
Under mechanical stresses and heat developed between screw and barrel during the extrusion process, the biopolymers undergo structural modifications: disruption of starch granules, melting of starch crystals, depolymerization of starch macromolecules, denaturation and aggregation of proteins, depolymerization of dietary fibers, and product browning due to Maillard reactions [3][4][5]. These changes are very important because they provide extruded starchy foods with an attractive appearance, structure, and relevant functional properties.
Due to an abrupt pressure drop at the die exit, the flash evaporation of superheated moisture forms a porous and expanded melt structure that is set after crossing its glass transition temperature. The structural features of the expanded foods govern their texture, functional and nutritional properties, and consumer acceptability [6][7][8][9][10][11][12]. The structural features consist of density, cellular structure, and properties of constitutive material, i.e., cell walls [13]. A phenomenological model of expansion allowed us to link extruded food features, such as expansion indices and cellular fineness to melt shear viscosity [14,15]. The melt viscosity takes into account the effect of formulation, such as fiber, protein, and starch contents, amylose to amylopectin ratio, and biopolymer structural changes during processing [16].
The development of mechanistic models of extrusion, based on the continuum mechanics approach, provides many benefits: a better understanding of the extrusion process and the development of computer-aided tools in simulation, optimization, and scale-up. In twin-screw extrusion modeling, two main approaches can be considered. In the first approach, local flow can be described with accuracy in a limited portion of the extruder by sophisticated numerical 2D or 3D models based on the finite element method (FEM) [17][18][19][20][21][22]. It can provide a very accurate description of the flow field only in a specific part of the extruder, and it is expensive in terms of computational resources and time.
The second approach is a global one; it is based on a simplified flow assumption and provides average values of the local flow variables (such as temperature, total dissipated energy, pressure, viscosity, filling ratio, and mean residence time) of the whole process, along the screws, from the hopper to the die exit [23,24]. It is based on various assumptions: (1) the flow is locally 1D or 2D; (2) the melting occurs instantaneously just before the first restrictive screw element; (3) the molten material behaves locally as a Newtonian fluid [24]. Computation of the flow parameters is done separately for each type of screw element and for the die components. The elementary models are linked together to obtain a global description of the flow along the extruder. This 1D global model has been implemented into a software called Ludovic ® v7.0.0 Classic Edition (Sciences Computers Consultants, Saint Etienne, France). It is quite straightforward to provide predictions of flow variables with simple computing resources and time: about a few seconds per simulation using a laptop configured with Intel ® Core™ i7-7600U processor (dual core). A fairly good agreement was obtained between simulation results from this 1D global model and a full 3D FEM model [18].
Despite significant progress in modeling and simulation, the design of extruded products at the industrial level is still based on a trial-and-error approach. One of the main challenges is determining the viscous behavior of melts under extrusion-like conditions that require specific rheometers, e.g., an in-line slit die rheometer or pre-shearing capillary rheometer (Rheoplast ® ) (Société Courbon, Saint-Etienne, France) [25,26]. This difficulty can be circumvented, in the first approach, by assuming that starch viscosity is a function of amylose content, and it can be interpolated from the published data on maize starches having several amylose contents [27].
The 1D global model has been widely used for the design of experiments of processing, but it is still limited to simple formulation (maize and wheat starches and wheat flour) [23,[28][29][30]. In process design, the model is exploited as a numerical screening tool of operating conditions to define and optimize the experimental domain in terms of flow variables (melt temperature, specific mechanical energy (SME), pressure). Its implementation in product design, for a wide range of complex formulations encompassing various structural and functional properties for one product, has never been tackled before. Until now, a simple link between a single property of extruded starchy products (for example, hydrosolubility) and a single process variable (for example, SME) was considered [29]. In this context, the aim of this paper is to test whether the 1D global extrusion model, implemented with an appropriate viscosity model, can be deployed as a prototype of a computer-aided tool for predicting properties of extruded foods and structural changes of biopolymers based on realistic formulations, from operating conditions. In other words, it will serve for product design applications.
For this purpose, we selected seven starchy food formulations to manufacture dense and expanded foods: cereal and tuberous-based pure starches and their blends with a wide range of amylose-amylopectin ratio, legumes-based ingredients (pea flour, fava bean flour, and fava bean starch concentrate), breakfast cereals, and blends of wheat flour-wheat bran. A detailed experimental database was built, including input and output data for extrusion simulation and extruded product structural and functional features. After describing the 1D extrusion model, virtual extrusion trials were performed, and validation steps were systematically detailed. The extruder operating curves representing a variety of predicted extrusion variables in the function of operating parameters were built. Then we strived to show how products with targeted structural and functional properties can be manufactured using the relations between extruder operating charts, predicted extrusion variables, and extruded product features, including biopolymer structural changes.

Building Database
Database for twin-screw extrusion simulation has been collected from literature and our own experimental results [31][32][33][34][35][36]. Due to the large number of experimental works, case studies were selected by applying several criteria: (1) presence of melt viscosity model and melting transition data as input, (2) existence of experimental data of melt temperature and specific mechanical energy (SME), and (3) properties of final products exhibiting clear relations with extrusion variables.

Raw Materials
Seven food recipes have been selected, reflecting a wide range of complexity of food formulations. The recipes were extruded to manufacture foods from the simplest ones, e.g., dense potato starch [32] and expanded starches (high amylose maize, waxy maize, and their blends revealing a broad variation of amylose content [27,33], and potato [31], to complex ones, e.g., breakfast cereals, and expanded snacks made from (1) legumes-based ingredients (pea flour [34], fava bean flour, and fava bean starch concentrate), and (2) blends of wheat flour and wheat bran [35][36][37] (Table 1). The extrusion data set of fava bean ingredients and breakfast cereals were produced from our original works. For these cases, raw material chemical composition was determined according to standard methods [34].
Two breakfast cereal foods were obtained by extrusion of a cacao-based formulation (called CB/Choco Ball) composed of a blend of rice flour, wheat flour, cacao powder, sugar, and salt and of a honey-based formulation (HB/Honey Ball) composed of a blend of corn flour, malt extract, sugar, and salt.
An overview of the chemical composition of all extrusion feeds is given in Table 1.

Extrusion Trials
The extrusion variables cover a broad spectrum of moisture content (15-36% wet basis), melt temperature (85-190 • C), and specific mechanical energy (100-5500 kJ/kg). An overview of extrusion systems for all feeds is given in Table 2, including equipment specificities.
Maize and potato-based expanded starches (Maize A, Maize B, Maize C, Maize D, and Potato 1; case n • 1) and wheat-based snacks (WF, WF-LB, WF-HB; case n • 7) were manufactured using pre-pilot extruders having a screw diameter between 25 and 55.5 mm [27,31,33,35,36]. The screw profile included conveying screw elements followed by a left-handed screw element. While maize starch blends were extruded through a long-slit die, the extrusion of potato starch and wheat flour-wheat bran used a short circular die. Pea flour extrusion (PF; case n • 5) [34] used a more complex screw profile: conveying elements followed by three kneading discs staggered at 45 • , and then three left-handed elements. Several die geometries were used, circular (D = 3 mm) of various lengths (5, 10, and 15 mm), attached or not with a multi-step slit die rheometer, as described in detail in Appendix A ( Figure A1).
Dense pea and potato starch (PS, Potato 2; case n • 3) [32] and fava bean-based snacks (FF, FSC; case n • 6) were extruded using laboratory-scale extruders with a screw diameter of 11 and 16 mm and equipped with a cylindrical die. The screw profile included conveying elements, followed by kneading discs staggered at 90 • and one left-handed screw element. For manufacturing dense (non-expanded) products (PS, Potato 2), the temperatures of the last barrel (Tb) and die (Td) were regulated at 80-95 • C, otherwise, the temperatures were set at melting temperature level.
Extrusion of breakfast cereal formulations (HB and CB; case n • 4) was carried out on a pre-pilot co-rotating twin-screw extruder (Coperion Werner & Pfleiderer ZSK 26Mc) equipped with a circular die. The same extruder and screw configuration, and temperature profile of barrel and die were used, as for pea snacks manufacturing [34]. The die had a diameter of 3 mm and a length of 15 mm. The extrusion simulation is based on a global model developed by Vergnes et al. [24]. This model uses a 1-D approach under the assumption that molten material behaves locally as a Newtonian fluid.
The flow is described by the Stokes equations in cylindrical coordinates (r, θ, z). Each conveying element is considered as a C-chamber associated with an angular sector representing the intermeshing area. Assuming that temperature and viscosity are locally constant, the velocity field can be determined in each element, and the volumetric flow rate Q v through the screw channel is expressed as: where A and B are coefficients depending on screw geometry, Ω is the screw rotation speed (rad·s −1 ), ∆P is the pressure variation (in Pa) along angular screw section ∆θ. The first term on the right side stands for drag flow, and the second one for pressure-driven flow. The non-Newtonian behavior of the molten material is taken into account by applying an appropriate model of shear viscosity η γ is the local shear rate. Specific models have also been developed for reverse screw elements [38] and blocks of kneading discs [39].
Flow equations are solved for each section of the screw and the die. These solutions are linked together to obtain a global description of the molten material flow along the extruder. The temperature changes ∆T in a section are obtained from the energy balance of accumulation, viscous dissipation E v , and heat transfer by conduction to the barrel E cd : where .
W is the power dissipated in the volume V of a screw element, ρ is the material melt density, and C p is its specific heat.
The melting zone is determined according to experimental observations: it is located at the first restrictive element (left-handed element, kneading disc), and melting is assumed to be instantaneous [40]. It is supposed that, at the point where the pressure starts to increase, upstream of the first restrictive element, the material is in the molten state and its temperature is equal to the melting temperature (T m ).
The computation starts from the die exit and proceeds backwards. As the final product temperature is unknown, an iterative procedure is used. A final temperature of the product at the die exit is chosen, and the pressure drop inside the die is calculated to go back to the upstream pressure value (to the upstream screw elements). The evolutions of pressure and temperature are thus determined element by element (section by section), from downstream to upstream until the melting point. The temperature calculated at this point is then compared with the melting temperature: if they are equal, the calculation is completed; otherwise, the final temperature is modified, and the calculation is repeated until convergence.

Input Variables
The input data are divided into three groups. The first group refers to extruder configuration. It covers screws, barrels, die, and also the position of specific zones (feeding, degassing). Examples of the screw profiles, for all feeds composed of right-handed and left-handed screw elements, and blocks of kneading discs are given in Appendix A ( Figure A2). Dies are characterized by elementary elements (e.g., slit, pipe, associated in series or parallel).
The second group refers to food properties: (1) Physical and thermal properties They cover density ρ, heat capacity C p , and thermal conductivity k in the solid and molten state, and melting temperature and enthalpy. Their dependence on moisture content (wet basis) can be calculated according to additive rules from the values for water and dry molten starches (Appendix A, Table A1) [41][42][43]. Most results on the melting temperature (T m ) as a function of moisture content of studied materials can be retrieved from the literature [32,34,42,[44][45][46][47] (Appendix A, Figure A3). The melting point of maize starches, having different amylose content (5%, 23%, 47%, and 70%), were obtained from the data set of waxy maize (amylose content of 5%); wheat (25%), smooth pea (35%), and wrinkled pea (75%) starches. The melting temperature dataset of breakfast cereal formulations (HB, CB), and fava bean flour and starch concentrate were determined experimentally as the peak of the DSC melting endotherm using standard procedure [32,48] (Appendix A, Figure A3).
(2) Rheological properties of molten starchy foods The variation of starchy melt viscosity (η) with shear rate . γ is generally described by the Ostwald-de Waele power law [25]: where n is the flow index (0 ≤ n ≤ 1) and K (Pa·s n ) the consistency.
The parameters K and n were expressed as a function of temperature (T, in • C), moisture content (MC), and SME according to [27]: n = n 0 + α 1 T + α 2 MC + α 3 SME + α 4 T MC + α 5 T SME + α 6 MC SME (6) where K 0 is the consistency at reference conditions, E is the activation energy, R is the gas constant, T a is the absolute temperature (in K), α is the water plasticization coefficient, β reflects the thermomechanical history coefficient, and n 0 and α i are constants accounting for the dependences of the flow index on T, MC, and SME and their interactions. The other constants, MC 0 , T 0 , and SME 0 , stand for reference moisture content (10%), temperature (353 K), and SME (350 kJ/kg). The viscosity model parameters of starchy melts, e.g., maize starches having different amylose content, potato starch, breakfast cereals, and blend of wheat flour-wheat bran were retrieved from the literature and are presented in Table 3. Under the hypothesis that viscous behavior is governed only by the amylose content of the starch, the viscosity of fava bean starch concentrate was approached by maize starch having an amylose content of 47%. The viscosity behavior of pea starch and fava bean flour was approached by the pea flour. Indeed, a previous study revealed that the master curve of shear viscosity of these feeds, determined by a pre-shearing capillary rheometer (Rheoplast ® ) (Société Courbon, Saint-Etienne, France), were very close to each other [49].
Finally, the third group of data concerns the operating conditions of the extruder. This data covers the barrel and die temperature profile, the moisture content, the total flow rate, the screw speed, and the heat transfer coefficients between molten material and barrel, screws and die, expressed through a Nusselt number. The data was overviewed in Table 2 for all recipes. For all simulations, regardless of the recipes, the Nusselt number was set to 20.
* The authors computed model parameters without taking into account reference variables (T 0 , MC 0 ). Therefore, the viscosity model took the form as;

Output Variables
The main variables predicted by the 1D global extrusion model (Ludovic ® v7.0.0 Classic Edition) (Sciences Computers Consultants, Saint Etienne, France) are melt temperature T, SME, pressure, and viscosity. These variables can be linked to the biopolymers transformation (e.g., starch melting and depolymerization, and protein aggregation), structural (e.g., expansion) and functional properties (e.g., water solubility index of starch) of extruded foods, considered as indirect output variables. Starch structural changes were represented by some features, such as the number of remnant starch granules/unit area, residual molar mass of starch, and intrinsic viscosity. Protein cross-linking by covalent bonds other than disulphide bonds was quantified by the amount of in-extractable proteins in sodium dodecyl sulphate (SDS) and dithioerythritol (DTE). Expansion behavior was characterized by density, sectional expansion index, and cell density (number of cells per volume unit). These properties have been listed in Table 4, hence completing the database, with their origin (literature, this work) and measurement methods. Actually, when these features were not available in the literature, they were characterized especially in this present work. In this case, methods have been detailed in Appendix A.2.

Validation of Extrusion Simulation
The reliability of the extrusion simulation was confirmed: (1) by comparing computed extrusion variables with the measured ones (T and SME) (

Comparison between Experimental and Simulation Results
Like in the experimental trials, the predicted melt temperature T_com was retrieved at the die exit. For most simulations, the computed temperature is overestimated by 5-10% compared to the measured temperature. The uncertainty of measurement can explain this difference. Temperature sensors are flush-mounted on the die wall with little penetration into the melt; consequently, the measured temperature differs from the actual one at the core of the product. Nevertheless, the computed temperature correlates well with the measured temperature, regardless of the formulation and extrusion system (R 2 = 0.88) (Figure 1a). Good correlations were also obtained between predicted and measured SME values (R 2 = 0.71-0.96), depending on the formulation (Figure 1b). The predicted SME is underestimated by 35-65% because the extrusion model only considers the melting energy and the viscous dissipation energy after the melting section, and neglects the solid transport and solid friction energy (inter-particle and particle/metal, metal/metal), which depend on raw material. All these trends in predictions: underestimation of SME and overestimation of temperature are in agreement with published results on extrusion simulation of starches and wheat flour [28,30]. The fair predictions of T and SME confirm the pertinence of the extrusion model for the design of extruded starchy products.

Axial Profile of Predicted Flow Variables
Simulation results were illustrated for two recipes of breakfast cereals taken as examples by presenting the axial profiles of the predicted variables (melt temperature, specific mechanical energy (SME_com), pressure, and viscosity) along the screws (Figure 2a) for typical moisture, feed rate and screw speed conditions (Figures 2 and 3). These profiles give us insight into thermomechanical conditions along the screws. The results obtained for other recipes exhibit similar trends and are not presented here.
The temperature increased progressively to values superior to the melting temperature T m (T m = 133 • C for HB and T m = 165 • C for CB) along the kneading discs and rose sharply along the left-handed elements (Figure 2b). The melt flow across the kneading discs caused a temperature increase of about 5 to 10 • C, depending on the material and flow conditions. A more severe temperature step was observed for left-handed elements: 35 to 60 • C for the three elements. The increase in temperature was essentially due to intense shear stresses across the restrictive elements causing high viscous dissipation. The temperature reached a maximum at the outlet of the last left-handed element. Then it decreased slightly in the melt conveying section before a final slight increase through the die.
The cumulative SME_com was calculated by summing values in each screw element. The melt flow across the restrictive elements (kneading discs and left-handed screw elements) increased the SME_com (Figure 2c), particularly in the left-handed elements. The sudden increase in SME_com was concomitant with the step increase in temperature, underlining the importance of viscous dissipation, which was mainly governed by the restrictive screw elements.
The axial pressure profile shows that the main part of the extruder was only partially filled since the relative pressure was zero for screw length between 50 and 200 mm (measured from the die) and beyond 300 mm. The left-handed screw elements and the right-handed screw element in front of the die were the only zones where the melt was under pressure (Figure 2d).
Since the viscosity depends on temperature and SME, the melt viscosity profile along the screws is also related to temperature and SME profiles (Figure 2b,c,e). At the melting section (in the first kneading bloc), the product temperature and SME were relatively low, and consequently, the local viscosity was high. In the following restrictive elements, and the conveying section nearby die inlet, there was an increase in temperature and SME due to the increase in viscous dissipation, leading to a decrease in viscosity. The knowledge of melt viscosity, mainly at the die exit, is very important, because it can be hardly measured experimentally and it partly controls the melt expansion at die exit [14].
The singularities of variable profiles (pressure peak, changes of temperature slope and energy increase and viscosity decrease) meet exactly the position of restrictive elements. This fact underlines the importance of viscous dissipation, mainly governed by the location of restrictive screw elements. Similar trends were observed for extrusion simulation of starches [28,29]. The CB recipe, once fully molten, exhibits higher temperature, SME, viscosity, and pressure along the whole screw length. Brugger et al. [50] observed that at close T and MC, the CB exhibited higher viscosity than the HB.
According to Figure 2b, the maximum temperature is predicted at the last restrictive element, i.e., at the position where it is difficult to measure it. For CB, at the selected operating condition (MC = 15%, N = 500 rpm, Q = 20 kg/h), the maximum temperature was high (230 • C) which suggests that food degradation may occur. For HB, under the same conditions, the maximum temperature was lower, but still very high for food stuffs (185 • C).
The importance of viscous dissipation on thermo-mechanical history was also proven by the negative effect of moisture content on the axial profile of all major flow variables (Figure 3a-c). When the moisture content is increased from 19% to 23%, the magnitude of temperature and SME_com is reduced all along the screws by a factor of 1.5 and 3, and pressure underwent a larger drop at the die entrance, by a factor of 4.5. The viscosity dependence on moisture content can explain this. Higher moisture content leads to a lower melt viscosity, which decreases energy dissipation and then lowers temperature and pressure.
Another important variable for a twin-screw extruder is the filling ratio, which is proportional to Q/N (Figure 3d). For each moisture content, the predicted SME_com of the CB recipe is negatively correlated with Q/N, following a power function. The same trend describing the effects of Q/N on SME was observed for experimental results of twin-screw extrusion of starchy products [28,32,34]. In fact, SME is usually proportional to N/Q: an increase of screw speed (N) or a decrease in feed rate (Q) will induce torque reduction, due to the shear-thinning behavior of melt, as well as the lower filling ratio of the extruder and more important viscous heating. According to Figure 3d, moisture content also affects SME-filling ratio curves negatively, which is explained by the influence of the water content on the viscosity.

Sensitivity Analysis
The impact of changing one operating conditions of the three ones (screw speed N, total flow rate Q, last barrel temperature Tb) on extrusion variables was examined systematically, while keeping other conditions (i.e., extruder configuration, screw profile) constant. Results are presented as 2D-operating charts and illustrated in Figure 4a-d for the manufacturing of breakfast cereal (HB recipe) at a moisture content of 19%.
As indicated previously, SME_com increases with screw speed and decreases with flow rate (Figure 4a). For each Q, SME increases by 20-30% when N increases from 150 to 700 rpm. At any N, SME drops by 35-40% when Q increases from 10 to 50 kg/h, with the highest decline occurring at the beginning of Q change (from 10 to 20 kg/h). The variations of T at die exit (T_com), when N and Q are changed, follow the same trend as SME_com, but it is quite insensitive to the change of Q (Figure 4b). The inverse trend with N and Q was observed for computed die entrance pressure (P_com): it decreases with N and increases with Q by a factor of 1.75-2.5 (Figure 4c). The negative impact of Q was also observed on predicted melt viscosity at the die exit: η_com drops by a factor of 2-2.25 ( Figure 4e). However, an increase of Tb from 90 to 140 • C leads to a decline in η_com by a factor of 1.2-1. 35. The increase of screw speed and the decrease of flow rate positively affects the viscous dissipation and, therefore, the SME and melt temperature. Indeed, the melt temperature increase leads to a decrease in the viscosity, and consequently, pressure.
Finally, the effect of conductive heating was studied on melt temperature. Regardless of N, an increase in barrel temperature (Tb), from 90 to 140 • C, leads to melt temperature (T_com) rising, by 5-10%, following linear trends (Figure 4e). However, the melt temperature is always higher than that of the barrel. It proves that viscous dissipation has more impact than the heat transfer between barrel and melt.
The set of Figure 4b-d can be used to evaluate the upper limit of operating pressure and melt temperature. For example, let us fix the upper limit of maximum melt temperature and pressure to be 165 • C and 10 MPa. The feasible operating domain for barrel temperature (Tb) of 90 • C at the highest flow rate (Q = 50 kg/h) is then screw speed in the interval of 300-550 rpm.
The order of magnitude of predicted extrusion variables and the trend of their responses to changes of extruder operating conditions and screw elements are similar to those observed for the extrusion of standard cereal and starches [28,29,56]. The similarity in the trend of responses was found for all studied starch-based recipes, likely because the responses are mainly governed by the viscous (shear thinning) behavior of the molten starch. These facts validate the 1D global extrusion model for the design of starchy foods. In addition, the 1D global model of twin-screw extrusion allows us to better understand the machine (extruder) and extrusion process itself. The knowledge of temperature profile and the possibility of adjustment of processing conditions using simulation appears thus advantageous, particularly for thermo-sensitive foods.

Application in Food Design
In this section, extrusion modeling is applied to determine the operating conditions that lead to a product with desired structure and properties. This approach is based on: (1) relationships between computed extrusion variables and product features (Figures 5a, 6a,b, 7a,b and 8a-c, and Appendix A, Table A2 for equations) and (2) using these relations and the global extrusion model to select an interval of operating conditions that permit to obtain the desired properties

Intrinsic Viscosity of Expanded Starches
Starch transformation was assessed by the intrinsic viscosity, reflecting the extruded product's molar mass. It was normalized to the property of native starch, resulting in a relative drop of intrinsic viscosity ∆Vis (%). The higher ∆Vis, the larger the extent of degradation, the more depolymerized extruded starch, and the lower the molar mass is. Extrusion conditions for expanded maize and potato starches (Maize A, Maize B, Maize C, Maize D, and Potato1) are indicated in Table 2.
For all recipes, an increase in computed SME (SME_com) led to a drop in intrinsic viscosity ∆Vis between 20 to 80%, which can be interpreted by starch depolymerization (Figure 5a). Conversely, no clear trend was observed for the impact of melt temperature (not shown). The trends are in accordance with experimental results (Appendix A, Figure A4a). The effect of viscous dissipation favors the disruption of starch granules and macromolecular chain splitting has been observed many times before, e.g., by Li et al. [3]. The relationships were manifested fairly well by power trends (R 2 ≈ 0.7-0.85), the parameters of which depend on starch origin (Figure 5a). Their numerical values are given in Appendix A, Table A2. It can also be seen that, in the SME_com interval [250, 500 kJ/kg], the values of ∆Vis rank inversely to their amylose content for maize starches. This is due to the larger molar mass of amylopectin, a highly short-branched macromolecule, which makes it more sensitive to shear, therefore to SME, than amylose, a linear macromolecule (10 8 vs. 2 × 10 6 g/mol). The extent of degradation ∆Vis is the most important in the case of potato starch, from 60% to 80% in the interval SME_com [250, 350 kJ/kg] (Figure 5a). This trend can be explained by the higher value of intrinsic viscosity for its native starch (450 mL/g) compared to the values given for normal maize starch (amylose content of 23%) in the native state (22.0 mL/g) [31]. Therefore, the larger size (molar mass) of potato starch macromolecules, amplified the sensitivity of this starch to shear.
Since it is possible to predict macromolecular degradation from a computed extrusion variable, it is interesting to predict the conditions for processing maize starch, e.g., having an amylose content of 47%, with a targeted intrinsic viscosity drop ∆Vis = 50 ± 5%. According to Figure 5a, SME_com should be in the interval [625, 875 kJ/kg]. Besides operating conditions (Q, N), viscous dissipation, reflected by SME, is sensitive to the profile and geometry of screw elements, mainly restrictive elements. Therefore, simulations were conducted to study the effect of the length of the left-handed screw element: (a) 0 mm, (b) 25 mm, and (c) 50 mm. According to the response surfaces, or operating chart (Figure 5b), at feed MC = 20%, when the restrictive element was absent (0 mm, case a), even by using most severe conditions (highest N and lowest Q), SME_com was lower than 310 kJ/kg, hence, lower than the target. In case b, when its length is 25 mm, the extruder can only operate in a very restrictive interval of Q and N values, represented by the hatched area: N ≥ 650 rpm and Q < 11 kg/h, which leads to SME_com lower than 665 kJ/kg, hence in the target interval (Figure 5c). Finally, for the longest restrictive element (50 mm, case c), the operating window was larger, with feed rate Q possibly reaching 20 kg/h for larger N values (750 rpm) (Figure 5d). This example underlines how the simulation can be used to select the relevant screw profile, which is helpful from a practical viewpoint, e.g., before implementing experimental trials.

Partially Melted Non-Expanded Starches: Study of Starch Structural Modification
Here, we deal with specific products, namely partially melted (partially gelatinized) dense extruded starches, i.e., Potato 2 and pea (PS). This kind of product was obtained through extrusion at the moderate moisture content (30%) and at the low barrel and die temperatures, i.e., 75-95 • C (Table 2, Appendix A, Figure A3). The starch transformation indicator was accessed at two levels: (1) granular one, expressed by the number of granular starch remnants per mm 2 (Ng), and (2) molecular one, revealed by the residual molar mass (RM, %) that is obtained by normalization to native starch property.
For both starches, the Ng, reflecting the extent of starch granule breakdown or starch crystallite melting, diminished with an increase in SME_com following exponential functions (Figure 6a). This phenomenon is common for starches extrusion [3]. In contrast with experimental results (Appendix A, Figure A4b), the Ng prediction of extruded pea starch was poor (R 2 = 0.56) and the fitting parameters depended on the raw material; this is likely due to the underestimation of SME prediction of pea starch extrusion. This inaccurate prediction comes from approaching the viscous behavior of pea starch melt with pea flour as an input of the 1D global extrusion model. Indeed, a previous study revealed that the shear viscosity of these feeds, determined by a pre-shearing capillary rheometer (Rheoplast ® ) (Société Courbon, Saint-Etienne, France), operated at SME < 250 kJ/kg, was very close [49]. Prediction improvement could be achieved by determining the viscous behavior of pea starch on a larger SME interval. For potato starch, the Ng dropped by 95%, while SME_com increased by a factor of 4 (R 2 = 0.73).
When the SME_com increases, a drop in RM occurs concomitantly with a decrease in Ng (Figure 6a,b). It means that the starch macromolecular chain splitting phenomenon is also involved in the starch crystal melting process; both phenomena being intensified with the more severity of shear stress [32]. The trends, specific for each product, can be fitted fairly well with power functions (R 2 = 0.7-0.89). The general tendency is in accordance with experimental results; however it is obviously seen experimentally that potato starch is more sensitive to molecular degradation than pea starch (Appendix A, Figure A4c). The RM of potato product dropped by a factor of 4 when measured SME was increased by a factor 1.5. On the contrary, the RM of pea product diminished by a factor of only 1.5 when the SME_mes was intensified by a factor 4. This trend was not visible for the relationship between RM and computed SME due to underestimation of SME prediction. The higher sensitivity of potato starch macromolecules to shear can be explained by (1) its higher amylopectin content, which leads to a higher molar mass and lesser tendency to entangle and (2) higher extrusion temperature: [T m , T m + 10 One example of a commercial application of this kind of product is dried noodles. The functional properties of noodles, such as optimal cooking time, water uptake, and cooking loss, depend on the state of starch transformation. During cooking (boiling), the remnant starch granules of noodles uptake water and swell, and cooking loss occurs through the diffusion (leaching) of depolymerized amylose and amylopectin outward the ghost starch granules, from the core to the surface of noodles. Cooking loss can be minimized by allowing starch swelling without excessive breakdown/leaching of its macromolecules. In other words, a compromise between two properties, Ng and RM, shall be sought, based on the assumption that the higher the RM, the more enhanced the leaching.
As an example, we want to manufacture partially gelatinized potato starch with the following features: (1) highly amount of undisrupted starch granule: 150 < Ng < 300; (2) moderate level of RM (RM = 50 ± 10%). According to the relationship described in Figure 6a,b, extrusion should be conducted at SME_com = 525 ± 25 kJ/kg. The operating chart of the extruder is given in Figure 6c. The feasible domain of screw speed and flow rate conditions is presented as a hatched area. The targeted product can be obtained at a low flow rate (0.6 kg/h) by tuning the screw speed N at 40-50 rpm. When a high throughput is required (Q = 1.8 kg/h), the screw speed should be set at a higher velocity: 140-165 rpm.

Legumes Based Snacks: Compromise between Protein Transformation and Expansion
Here, the features of protein-rich snacks, extruded from pea and fava bean ingredients, were represented by protein solubility, in terms of insoluble proteins and density. While the latter reveals the macrostructure, the first can reflect the structure at supramolecular level. The gain in insoluble proteins indicates a higher amount of proteins crosslinked by covalent bonds other than S-S ones, i.e., isopeptide bonds. The higher is the density and the less expanded is the snack.
The rising of melt temperature, at die exit, observed in the pea flour extrusion simulation promoted the protein cross-linking by non S-S bond (R 2 ≈ 80, Figure 7a), in accord with experimental results (Appendix A, Figure A4d). The onset of protein cross-linking occurred at a computed melt temperature (T_com) of 155 • C. The snack density decreased from 800 to 100 kg.m 3 , in other words, expansion was magnified by a factor of eight, when T_com increased from 125 to 175 • C, regardless of formulation (R 2 = 0.75, Figure 6b). The positive impact of temperature on legumes expansion, in agreement with experimental results (Appendix A, Figure A4e), is not common for extrusion of starches, where the inverse trend is usually observed [14]. The expansion mechanisms are likely governed by melt rheological changes linked to a certain morphology of starch-protein composite that is probably tuned by structural protein modification [34]. Numerical values of all fitting parameters are summarized in Appendix A, Table A2. The impact of SME_com on snack density and protein transformation is not clear (R 2 = 0.45 and 0.55, result not shown).  A low density (lower than 200 kg/m 3 ), resulting from high-temperature extrusion (155 < T_com • C < 175 • C) at low MC (18% and 21%), is indispensable for obtaining crispy snacks. However, extrusion at these temperatures will be accompanied by a higher amount of insoluble proteins (from 5 to 25%). In addition, according to some authors, the protein cross-linking can deteriorate nutritional properties, such as, i.e., protein digestibility [10], and cause a decline in customer acceptance because protein cross-linking increases viscosity of chewed food [9], leading to more difficulty in food swallowing. Therefore, the design of protein-based snacks like pea snacks requires a compromise between density and protein cross-linking.
An example of compromise is set as follows: snacks having a density lower than 200 kg/m 3 and insoluble proteins lower than 10% can be extruded at computed temperature (T_com): 155 < T_com • C < 165 (see correlation charts, Figure 7a,b). The operating charts of extruder for extrusion at MC of 21% are given in Figure 7c

Expansion of Wheat Snacks Enriched with Wheat Bran (Fiber)
The features of wheat-based snacks are represented by hydro-solubility of starch (WSI starch ), macrostructure described by sectional (radial) expansion index (SEI), and cellular structure expressed as cell density per cm 3 (NC). The higher the WSI starch , the more soluble are the starch molecules. This variable can thus indicate the state of starch depolymerization induced by extrusion. The higher the NC, the more porous the expanded structure.
Regardless of fiber (bran) content, all wheat flour-based snacks exhibited a unique linear and positive relationship between WSI of starch and SME_com (Figure 8a), in agreement with experimental results (Appendix A, Figure A4f). The higher the shear energy, the higher the extent of starch depolymerization, which is a common trend in starch extrusion [3]. More soluble and smaller starch molecules are more accessible to digestive enzymes during food digestion; however, a highly digested starch can lead to a high value of the glycaemic index.
The prediction of WSI starch from both computed and measured SMEs is poor (Figures 8a and A4f). If the fiber (wheat bran) is water-soluble, the data dispersion can be likely due to competition in water absorption between starch and fiber; otherwise, the dispersion origin can be explained by filler content (insoluble fiber). The increase of filler content in continuous molten starch can increase melt viscosity, affect viscous dissipation energy, and finally, the degree of starch transformation in the extruder [16]. The prediction could be improved by several means: First, presenting the WSI of total solids, considering fiber and starch, instead of the WSI of starch only. Second, determining the quantity of soluble starch experimentally using a more sophisticated method, such as orcinol titration [34] rather than the gravimetric-based method used by Robin et al. [36]. Third, improving the viscosity model, the effect of filler volume fraction could be taken into account.  The melt temperature and viscosity at the die exit was predicted by simulation since experimental data was not available. Acceptable correlations following exponential functions (R 2 = 0.7-0.8) were found between T_com and structure, such as SEI (indicator of expansion in the radial direction) (Figure 8b) and cell density (porosity indicator) (Figure 8c). However, the bran inclusion led to a similar trend but different correlations compared to the flour without bran. We observed that melt temperature and bran enrichment negatively affected SEI. In contrast to SEI, bran inclusion and melt temperature positively affected cell density (NC). For bran-enriched snacks, NC escalated by a factor of 10 when T_com was raised from 125 to 195 • C. The onset of porous structure formation took place at close T_com (155 and 160 • C) for snacks made without and with bran inclusion, respectively. At this onset temperature, bran inclusion surged NC by a factor of 6. Wheat bran can be considered as a filler that can enhance nucleation and hence expansion [16]. An increase in melt temperature intensifies bubble nucleation and bubble growth and thus promotes expansion and porous structure (high NC) [14]. The cell density also negatively correlated to predicted melt viscosity (η_com) (R 2 > 0.95, Figure 8d), with distinct power function expressed to each bran content. This trend can mainly be attributed to the negative effect of shear viscosity on the expansion mechanisms such as nucleation, bubble growth, and setting [14]. Higher porosity can make a snack texture more fragile/crispy. Therefore, the design of wheat-based snacks requires a compromise between hydrosolubility of starch and porosity of expanded structure.
In general observation, the 1D global extrusion model is more accurate in the prediction of extruded legumes features than extruded starch ones. The protein structural changes of extruded legumes were predicted from the computed temperature at die exit (Figure 7a,b), while the starch structural changes from computed SME (Figures 5a, 6a,b and 8a). The overestimation of temperature prediction by the extrusion model is much lower than the underestimation of SME prediction (5-10% vs. 35-65%) (Figure 1) and, therefore, a better prediction of extruded legumes features can be expected. The extrusion model neglects the solid transport and solid friction energy and material sticking, leading to SME underestimation. The friction and sticking phenomena depend on the raw materials. It explains why the prediction of intrinsic viscosity of extruded potato starch, through computed SME, is poorer than that of other starches (Figure 5a). Better SME prediction, and thus better product features prediction, could be achieved by improving the extrusion model by taking into account the friction and sticking coefficient.

Conclusions
Twin-screw extrusion processing of starchy products was simulated using the 1D global extrusion model (Ludovic ® software v7.0.0 Classic Edition (Sciences Computers Consultants, Saint Etienne, France)), and melt rheological model as input. A large data set has been built, including results for extruded products' properties and extruder variables (temperature and specific mechanical energy), taken as output variables. Satisfactory correlations between process parameters predicted extrusion variables and features of extruded foods, such as structural and functional ones, justify the use of the extrusion model as a computer-aided tool for designing starchy food by extrusion. Moreover, the capacity for predicting these features is extremely vital to reducing time and labor costs in industrial foods R&D.
Flow variables at the die exit, such as melt temperature, specific mechanical energy and melt viscosity, can be computed and used, in turn, as the input variables for a phenomenological model of expansion, allowing predicting expansion indices and cellular structure of expanded starchy foods. The possibility to compute melt viscosity is particularly interesting since it is difficult to measure, either in-line or off-line and because it governs the expansion phenomenon, at least for low hydrated starchy foods. In addition, the knowledge of temperature profile and processing conditions resulting from the simulation is useful in setting extruder configuration and operating parameters for extrusion of thermo-sensitive foods.
As a future prospect, extrusion simulation will be applied to the design of a wide range of starch-protein blends and expanded snacks from pulse crops. Further possible extension includes the processing of plant proteins for meat analogues. Besides classical structural and functional properties, the textural and nutritional features, such as starch and protein digestibility, inactivation of anti-nutritional factors, and destruction of pathogenic microbes and vitamins, will be tackled. The rheology database will be updated with melt viscosity models encompassing a wide range of moisture content, temperatures, and SME, taking into account protein and fiber content, as well as protein structural changes and starch-protein morphology.

Data Availability Statement:
The data presented in this study are available in Appendix A, and the dataset can also be extracted from the Results Section (Figures 1-8, Tables 1-4).

Acknowledgments:
We are thankful to Imen Jebalia (INRAE) for the density measurement of fava bean based-snacks.

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1. Figure A1. Die profile used in pea flour extrusion. It consists of a combination of a rheometric die and a cylindrical die. The schema is drawn with Ludovic ® software according to description given by Horvat et al. [26]. The rheometric die is a slit-die built up by a triple-step insertion module in a slit channel having a width of 30 mm. The geometry (length and height) of each step is given as follows: step 1-170 × 7 mm 2 ; step 2-120 × 5 mm 2 , and step 3-80 × 3 mm 2 . The cylindrical die is built up by a plate with a hollow/cylindrical channel having a diameter of 3 mm and three options of length (5, 10, and 15 mm). (d) Figure A2. Figure A2. Illustration of extruder screw profiles for manufacturing of (a) Expanded sample from maize starches having different amylose content [27,33] and potato starch [31]. For potato starch, two cylindrical dies in parallel were used, (b) Dense sample from potato and pea starches [32] and expanded samples from fava bean flour and starch concentrate, and from blends of tapioca, Empure and Merizet starches, (c) expanded snacks from pea flour [34] and breakfast cereals products. For pea flour, the presented extruder was equipped with a cylindrical die only. (d) Expanded snacks from bran-enriched wheat flour [35,36]. All screw representations were build-up using Ludovic ® software according to description given in the literatures.  Figure A3. Inverse of melting temperature (T m ) as a function of water volume fraction v. The points refer to melting temperature obtained as the maximum temperature of an endothermic peak of DSC. The curves represent the best fit using the Flory-Huggins equation [57].