A Review of Dynamic Models of Hot-Melt Extrusion

Hot-melt extrusion is commonly applied for forming products, ranging from metals to plastics, rubber and clay composites. It is also increasingly used for the production of pharmaceuticals, such as granules, pellets and tablets. In this context, mathematical modeling plays an important role to determine the best process operating conditions, but also to possibly develop software sensors or controllers. The early models were essentially black-box and relied on the measurement of the residence time distribution. Current models involve mass, energy and momentum balances and consists of (partial) differential equations. This paper presents a literature review of a range of existing models. A common case study is considered to illustrate the predictive capability of the main candidate models, programmed in a simulation environment (e.g., MATLAB). Finally, a comprehensive distributed parameter model capturing the main phenomena is proposed.


Introduction
Hot-melt extrusion has been a forming technique in popular use in industry since the 19th century.It is a thermomechanical manufacturing process where solid materials are transformed into a uniform product with a specific, possibly complex, shape by forcing passage through a die [1].Besides the die, the process includes a conveying system, involving the rotation of a screw (or several screws), such that the material is transported by an Archimedean screw action through the device.Nowadays, more and more forming processes use this technique, and the products that can be extruded are numerous: plastics [2,3], drugs [4][5][6][7], food [8][9][10], metals [11,12], tissue engineering scaffolds using electrospinning [13][14][15], porous items using supercritical fluid [16], etc.
A discontinuous extrusion process, which is also widely used in industry, is the ram extruder, where the whole material billet is directly injected into the device and pushed out from the die by the ram pressure [17].This single or multi-ram batch process can generate extremely high pressures, but usually results in poor melting efficiency [18].In the following, this type of process will not be discussed, and attention will be focused only on modeling twin screw systems.
Hot-melt extrusion is a complex process, as it involves an assembly of different screw elements (forward element, reverse element, kneading block) to manufacture the end product.The selection of the design parameters of the screw and the operating conditions of the process, including the screw rotation speed N, the barrel temperature T b , the feed flow rate of material Q in and the component concentration C in (see Figure 1), offer a wide range of possibilities.In turn, this multivariable process is characterized by a set of internal state and output variables, including material filling ratio f , pressure P, material temperature T m , screw temperature T s , die flow Q out , etc. (see Figure 1).To regulate the process around the desired operating point, a control system is usually required [19][20][21][22][23][24].The aim of this paper is to review the mathematical models of hot-melt extrusion that were proposed in the past few decades, with a focus on dynamic models that could be the basis for developing a simulator, a state estimator or a controller.The reviewed models range from simple continuous-stirred tank reactors (CSTRs) in series to sophisticated distributed parameter models, involving mass and energy partial differential equations (PDEs).The interesting approaches of population balance modeling (PBM) and discrete element modeling (DEM) are addressed, and purely data-driven, black-box models, are also considered.
In order to not restrict the review to a purely theoretical descriptive point of view, a simple case-study is considered to benchmark the main model candidates.To this end, dynamic simulators are developed in the MATLAB environment.As a result, a dynamic PDE model is proposed as a relevant model of the extrusion process.
This paper is organized as follows.The next section introduces the case study that will be used to illustrate and test the main modeling approaches.Section 3 presents several explicit models describing the residence time distribution (RTD) in the extrusion device, while Section 4 discusses dynamic models based on a discretization of the system in CSTRs and differential equation (ODE) models.Partial differential equation (PDE) models are then reviewed in Section 5. Section 6 provides a brief introduction to PBM and DEM, while a review of black-box models, i.e., transfer function models and neural network models, is presented in Section 7. The last section of this paper draws the main conclusions.

A Simple Case Study
A simple case study is used as a benchmark in this article, so as to illustrate the main modeling paradigms.It consists of a twin screw with the same pitch all along the device (see Figure 2).The geometrical and physical parameters are listed in Table 1.Moreover, this table gives the parameters of the empirical equations selected to model the material rheological behavior.Various studies were dedicated to the evolution of material characteristics, such as the viscosity, melt density, size and shape distribution of the particle clusters, electrical properties and evidenced complex dependencies with respect to the degree of mixing, the quantity of air into the extruded compound, etc. [25][26][27][28][29][30][31][32][33].For some extruded materials, analytical studies do not allow one to obtain an accurate representation of the rheological properties, and in-line measurements from an actual extrusion device are required [34].Rheological properties largely contribute to the nonlinear process behavior.In our case study, only the shear viscosity variations are considered, and Yasuda-Carreau equations, which are functions of the shear rate and the temperature with four adjustable parameters, are employed (see, for instance, [35][36][37]).They are described by the following equations: where η 0 (t) is the zero shear viscosity, λ t the characteristic matrix time constant, γ the shear rate, a the Yasuda parameter and n p the pseudo-plastic index.The dynamic zero shear viscosity η 0 (t) is varying with the material temperature T following the Arrhenius law:  Using one of the most elaborate models [38], which is presented in detail in Section 5, a simulation code is implemented in MATLAB.The partial differential equation model is discretized using four-point finite differences with 150 nodes, and time integration is performed with the solver ode15s.A residence time distribution (RTD) test is simulated with a screw speed rotation N = 100 RPM (revolution per minute) and a feed flow rate Q in = 0.358 kg/h.An additive Gaussian white noise with a standard deviation of 0.002 is considered, to represent the situation where the RTD is measured using tracer and colorimetric techniques [39].The resulting RTD is represented in Figure 3.

The Earlier Models: Static Maps
Basically, extrusion enhances component mixing.A fundamental process characteristic is the time spent in the extruder by the several components of interest, which can be assessed thanks to an RTD, as shown in Figure 3.According to [9,40], the RTD signal reflects the process operating conditions, such as the material homogeneity, the uniformity degree of particle shearing, the calorific exchanges (input barrel temperature or outside thermal flux changes) depending on the screw rotation speed, the feed flow rate, etc.The capacity of a model to predict the RTD of each component is therefore important to check that mixing, and possibly the chemical reaction, is achieved before the die.The analysis of the RTD shape leads to a general description of the flows inside the device.Consequently, monitoring this output is particularly important during process scale-up [41].
Earlier modeling works focused on an explicit representation of the RTD (see [42,43]).The first extrusion models, as reviewed in [44], consider and combine two types of reactors with a priori known RTDs (see Figure 4): • Plug flow reactors characterized by a uniform residence time distribution for all particles.
• Continuous stirred-tank reactors (CSTRs) assuming the non-uniformity of the residence time, but a uniform composition of the flow.Based on these building blocks, the first model consists of a cascade of CSTRs (see Figure 5), while the second model structure combines a plug flow reactor with a cascade of CSTRs.The mathematical model corresponding to these structures can be expressed in the form of Equation (3): where p = t min t is the plug flow effect and b = n (1−p) is the ratio between the diffusive and the convective effects.t min is the minimum residence time; t is the mean residence time; n is the number of reactors; and τ = t t is a dimensionless time variable.Parameter identification based on RTD experimental data is straightforward, but the model fit is not always satisfactory [45], particularly with the plain CTRs-in-series model.The plug flow model compartment is often required and can improve the data reproduction to some extent [45,46].
The selected number n of reactors varies significantly with the application, from eight in chocolate manufacturing [46] to 500-900 for a granulation process [45].The reader is referred to the works [47,48] where the appropriate number of reactors is determined for different screw elements based on a thermo-mechanical partial differential equation (PDE) modeling.
These preliminary observations demonstrate that the physical phenomena occurring in an extrusion process are difficult to model.For instance, the die shape causes obstructions, resulting in leakage flows and material stagnancy.Efforts were therefore deployed in order to consider these aspects (mathematical expressions are not included in the text for the sake of conciseness): • The work in [49] models the RTD of wheat flour particles using five volume elements of different sizes, each sub-divided in a series of CSTRs with identical volumes.In order to represent the tail in the RTD, three reactors with backward flows are included so as to explain material accumulation at the die (see Figure 6).• The work in [50] uses a combination of a plug flow reactor with a cascade of CSTRs with the same volume and a backward flow Q B (see Figure 7).The plug flow reactor represents the feeding zone and is characterized by a delay t d .The cascade of CSTRs models the completely filled part inducing a reflux and characterized by two parameters: the number of reactors n and the ratio between the forward and backward flows σ = Q Q B .• The works in [51] and [45] suggest another model inspired by [44] for a granulation process.
A plug-flow reactor is placed in series with a cascade of CSTRs with dead zones (see Figure 8) where a fraction of the material is retained.The earlier models (which at least find their roots in earlier works, some of them having been proposed relatively recently) are therefore static maps based on RTD calculation reflecting the operating conditions (screw configuration, screw rotation speed, feed flow rate, etc.) linked to the desired end-product quality.
To illustrate their potential, the model of [51] is used to represent the RTD measurements of Figure 3. Equation ( 3) is modified to include a dead zone parameter d, which influences parameter b = 1 A least-squares cost function is defined in order to estimate the unknown parameters θ = [n, t min,d ] through the solution of the following problem: The estimated parameters are listed in Table 2, and a model direct validation is presented in Figure 9, which is quite satisfactory.Static maps are limited in scope and do not offer much potential in terms of process control.The next section therefore pays attention to differential equation models, based on similar decomposition of the extrusion device into specific elements, such as CSTRs.

Reactors-In-Series Approach
The reactors-in-series approach is a direct extension of the previous concepts, i.e., it is based on a decomposition of the extrusion device into volume elements where perfect mixing is assumed, but dynamics are now taken into account.Based on mass and energy balances and on rheological laws, variations of the main state variables, such as filling ratio (ratio between the volume occupied by the material over the available mixing volume; a value of one represents a completely filled screw channel section and less than one a partially filled screw channel section), pressure, temperature and monomer concentration, can be represented by ordinary differential equation systems (ODEs).
Mass balances essentially allow one to determine the filling ratio of each reactor based on the inlet (Q in ) and outlet (Q out ) mass flows.Considering n reactors, mass balances are written as follows: where i is the reactor index, ρ is the material density and V i and f i are the respective volume and filling ratios.
Mass balances can also be used to determine the reacting monomer concentration, as in: where M is the molar mass of the monomer, C in the inlet monomer concentration, C i the concentration of monomer in reactor i and R the reaction rate.
Temperatures can be determined on the basis of energy balances involving several heat transfer mechanisms: convective heat transfer, heat exchange between the screw elements and chemical reaction heat.In many works [35,[52][53][54][55], viscous friction dissipation is also included in the energy balance.The computation of the barrel and screw temperatures is possible, but for the sake of clarity, only the expression of the material temperature T m is given here: where M m i is the mass of material in reactor i, c m p the specific heat capacity, U b and U s the heat exchange coefficients between material/barrel and material/screws, respectively, S b and S s the exchange surfaces between material/barrel and material/screws, respectively, r the reaction rate and ∆H the reaction enthalpy.
Since density and shear viscosity enter into the expressions of the flow and viscous friction dissipation, rheological properties of the material should be included in the process model.In the literature, a constant density is often used and the shear viscosity is modeled by an empirical equation [56][57][58] or a rheological law [35,[52][53][54][55].
The reactors-in-series models can be classified into two categories depending on whether specific sections of the screws (C-chamber, inter-meshing element, kneading block) are represented by a CSTR.A classifying review of the existing models is achieved in the following.

• Models with Predetermined Screw Discretization
In [56,57], the filling ratio, monomer concentration, temperatures and RTD are determined for a counter-rotating twin screw extruder using a subdivision of the screws in C-chambers (see Figure 10) where mass and energy balances are expressed using continuum mechanics.If the system includes a high inter-meshing area between the screws, other volume elements may be added to model the specific flows in this zone (see [52]).An extrusion device presents a varying filling level, a function of the local screw geometry.It is generally considered that an extruder comprises two different zones, which can be partially or completely filled, inducing specific mass balance expressions.Indeed, mass flow networks are different in each zone.In partially-filled zones, only a pumping flow Q p , due to the screw rotation, transports the material along the section.In completely filled zones, the presence of a pressure gradient creates leakage flows, and the resulting mass flow network is modeled by parallel C-chambers, in which the leakage locations are divided into five groups [59] (see Figure 11  In these models, pressure can be considered as a system input in completely filled zones [56,57] or can be determined by integration of the filling ratio, as in [52].

• Models with Adjustable Discretization
In [35,53], the different screw sections can be discretized using an adjustable number of CSTRs (see Figure 12).Material can flow in the forward direction, with flow Q f due to the screw rotation, and in the backward direction, with flow Q b due to the pressure difference between two reactors, as shown in Figure 12 and Table 3 [60].Table 3. Expressions of the forward and backward flow rates where K f and K b are screw geometric parameters, P is the pressure, f the filling ratio and V the reactor volume.
Based on the filling ratio, the pressure in each reactor can be computed using a simple algebraic method.Two different situations may happen: in partially-filled reactors ( f < 1), the pressure P matches the atmospheric pressure P 0 , and in the completely filled reactors ( f = 1), the pressure is above the atmospheric pressure; and flow continuity is verified.The problem can therefore be expressed as a system of linear algebraic equations Ax = b, where x = P is the pressure vector and A a tridiagonal matrix.
Some additional improvements are proposed in recent models, such as refined expressions of the kneading block flow in [54] or a new geometric definition of parameters K f and K d in [55], as well as the assumption that flows inside the kneading block are only described by pressure gradients.
As an application example, the model proposed in [35] is selected to represent the data provided in Figure 3. Three parameters have to be estimated: the number of reactors n and two flow geometric correction factors A f , A b , for the forward and backward flow (these latter parameters are multiplicative factors of the screw geometric parameters K f and K b and take into account the assumption that screw channels are rectangular in shape).
The estimated parameter values are given in Table 4, and direct validation is illustrated in Figure 13, which is obviously a better fit to the experimental data than Figure 9.An additional advantage of ODE models is the possibility to simulate various scenarios.For example, the effect of step changes in the inputs can be easily studied, e.g., a step change in the screw speed rotation N from 100 down to 75 RPM and a step-change in the feed flow Q in from 0.358 down to 0.15 kg/h.The choice of Q in as an input includes the use of the "starve-fed" mode where the flow is precisely controlled by loss-in-weight feeders.It should be noted that "flood feeding", where the screw rotation speed controls the flow rate, can be another method to produce feeding profiles.The process outputs, e.g., the die pressure P end and the outlet flow Q end , are presented in Figures 14  and 15, respectively.These results show that dynamical analysis of the extrusion system can be easily performed.However, model accuracy is still perfectible in several aspects.For instance, poor agreements between simulation and experimental data at a low screw rotation speed have been reported in [56] or a discrepancy of the influence of the screw rotation speed on the pressure die has been observed in [52].This motivates the exploration of distributed parameter models, involving mass and energy balance partial differential equations.

Distributed Parameter Models
The direct formulation of mass, energy and, possibly, momentum balance partial differential equations (PDEs) is an appealing alternative to the use of the CSTR-in-series paradigm, where there is an interplay between the model definition and the numerical solution via the number of reactors.Distributed parameter models allow the uncoupling between model formulation and solution techniques, which can be of arbitrary sophistication (finite differences, spectral methods, finite volume methods, etc).Depending on the modeling assumptions, i.e., 1, 2 or 3 spatial dimensions, a (non-)Newtonian fluid, (non-)isothermal conditions, many different model formulations are possible.
3D models are presented, for instance, in .These models are solved using classical computational fluid dynamics (CFD) tools, such as finite element or finite volume methods.
These models can include all of the features already discussed in this article, such as a description of the changes in the rheological properties of the material along the device, the presence of partially-or completely-filled zones, the wall slip effect on the flows [83][84][85][86][87], etc.The resulting 3D models can be used to achieve some dynamic studies, such as the analysis of the instabilities of the die flow [88,89] or the description of the chaotic motion inside screw elements by dynamic tools (Lyapunov exponents, Poincare sections and particle tracking), which enable the study of the mixing conditions in hot-melt extrusion processes [90][91][92].However, the numerical solution of the conservation laws requires important computational effort, and 3D CFD will be mostly useful for system design and steady-state analysis.
The work in [93] presents a comparison between a 3D and 1D PDE model of a specific extrusion device.For dynamical study and model-based control, 1D representations will usually be more appropriate.Among the works considering a 1D axial representation, one can note [94][95][96][97][98][99][100][101][102].In these works, the spatial profiles of the filling ratio, pressure and temperature along the different screw elements are studied in steady state.However, most of them could be applied for transient computation, as well.
The literature shows that many 1D dynamical models are based on [38], in which the system is divided into two specific zones: the solid conveying zone (SCZ) and the melt zone (MZ) (see Figure 16).Steady state analyses of such models are achieved in [103].-Mass balance: where ξ and N are respectively the screw pitch and rotation speed.-Energy balance: where T is the material temperature, φ the viscous energy dissipation factor, c s a geometric parameter characterizing the viscous energy dissipation [104], η s the viscosity, ρ the density, V the shear volume, c p the specific heat capacity, U s the convective heat transfer coefficient, c e the equivalent circumference of the twin screw and T b the barrel temperature.
• Melt zone: Unlike the solid conveying zone, the forward flow is considered as constant since the density itself is constant and the zone completely filled.In the following equations, the backward flow is associated with the pressure flow between two parallel plates assuming that the screw gaps are much smaller than the screw diameter.The net flow rate at any location in the melt zone equals the output flow at the die, Q d .The mass and energy balance equations are: -Mass balance: where B is a constant governing the leakage back flow and P the pressure.
Determination of this pressure is possible using the following equation: -Energy balance: where c m is a geometric parameter characterizing the viscous energy dissipation [104] and η is the viscosity.U is the convective heat transfer coefficient.
Due to the existence of two zones, a boundary condition at the moving interface must be expressed.During the transient period, the melt zone length (l m ) changes.This can be explained by an overall mass balance of the melt zone: where l m is the melt zone length, f s is the filling ratio at the melt zone inlet, Q s is the inlet mass flow rate in the melt zone and Q d is the output die flow rate.2V(1 − f s ) represents volume free from material.An ODE is formulated to determine the melt zone length: This dynamic model results in a PDE system coupled to an ODE describing the moving boundary between the solid conveying zone and the melt zone.The movement of this separation involves the zone description (solid conveying or melt), its discretization and the associated balance change during time.The application of standard numerical techniques is therefore difficult, and an alternative method is preferred.The principle of a simplified solution procedure is represented in Figure 17.The suitability of this model [38] has however not been tested in a more realistic screw configuration (e.g., more complex geometry, presence of kneading blocks, etc).In real applications, the screw geometry discontinuity involves a discontinuous filling ratio f evolution, which can cause numerical problems.The authors of [105] alleviate this by focusing on a continuous variable: the net flow rate Q.In a configuration with reverse pitch sections, they propose a mass balance equation allowing the determination of Q: To tackle the problem of a moving interface between the partially-and completely-filled zones, [106] suggests the introduction of a new spatial variable χ ranging from 0-2 (see Figure 18).The solid conveying zone is delimited by the interval [0,1], and the spatial variable is expressed as χ(x, t) = x l s (t) .In the melt zone, the new variable is restricted to the interval [1,2] by χ(x, t) = x+L−2ls(t) L−l s (t) .These changes of variable are convenient as they transform a PDE problem with a moving boundary into a PDE problem on fixed spatial domains, which makes spatial discretization easier to implement.In this review, a dynamic model is implemented based on these recent works and, notably, the work in [38].The filling ratio, pressure and temperature are modeled using a division of the device in two zones, and an ODE is added to the equation system to determine the moving interface.Using the method of lines (MOL) [107] with a four-point finite difference scheme on a spatial grid with 150 nodes and the ordinary differential equation solver ode15s, numerical solutions can be computed accurately.To reproduce a tracer concentration, such as the one used in the RTD test, Equation ( 16) is used: where C is the tracer concentration, Q the net flow rate and D the diffusion coefficient.The data provided in Figure 3 is used to estimate the parameters in Equations 8-14 and 16, e.g., the shear volume V, a constant governing the leakage back flow B and the diffusion coefficient D. These unknown parameters are identified using a least-square cost function similar to 4. Estimated parameter values are summarized in Table 5, and a direct validation is presented in Figure 19.The estimated parameters actually correspond to those used in the original simulator, which generated the RTD "experimental" data.Confidence intervals shown in Table 5 demonstrate that the information content of the RTD signal is sufficient to ensure the practical identifiability of the three considered parameters.As with the CSTR-in-series model, the time responses to step changes in the input can easily be evaluated, i.e., die pressure P end and the outlet flow Q out responses to step changes in the screw rotation speed rotation N from 100 down to 75 RPM and in the feed flow rate Q in from 0.358 down to 0.15 kg/h.The results are shown in Figures 20 and 21   This dynamic analysis leads to output trajectories comparable to the previous results (see Figures 14 and 15).However, a slower transient response is observed for a feed flow Q in step in the case of the PDE model.

Population Balance Modeling and Discrete Element Modeling
This article mostly focuses on dynamic modeling of the global characteristics of the twin screw extruder (filling ratio, pressure, temperature).However, some manufacturing processes require specific material properties for their implementation, and the knowledge of their dynamic evolution along the extruder is important.Population balance modeling (PBM) and discrete element modeling (DEM) are two methods dedicated to this kind of problem, which are briefly discussed in this section.
In PBM, the individuals (particles) are characterized by several properties (independent variables), such as size, liquid content, porosity, spatial coordinates, etc, and are involved in several rate processes, such as aggregation, breakage, nucleation and growth.The models consist of integro-differential equations, in one or more dimensions, which can be formulated as [108]: where F(x, t) represents the number of particles as a function of particle class x and time t.R f ormation (x, t) and R depletion (x, t) are respectively the rate of formation or depletion of the particles depending on several rate processes.These rate processes are generally represented by empirical or semiempirical expressions (for instance, [109][110][111][112]), which is one of the main disadvantages of this modeling approach.Experimental data are required to develop and validate these expressions in some restricted operating regions.PBM are used in various areas, such as crystallization [113], powder mixing [114], powder milling [115] and wet granulation [116,117].In extrusion processes, this method can be employed to compute different material properties, such as the active pharmaceutical ingredient (API) or excipient component concentration [118], the RTD [119], the particle size distribution, the liquid content and the porosity [120,121].
DEM is based on a different approach where each particle is considered as an individual entity [122][123][124][125][126][127].Two classes can be defined: hard-sphere methods where particles are considered rigid and collisions are instantaneous and binary (note that this method is not valid in twin screw modeling) and soft-sphere methods where interactions upon collisions and lasting contact between the particles and the equipments (screws, barrel) are considered.The implementation of the collisions and contact forces (normal and tangential forces; see Figure 22) in Newton's second law of motion (Equation ( 18)) and the angular momentum (Equation ( 19)) equations at any time instant constitute the modeling.
where v is the particle speed, m the particle mass, g the gravity, I the particle inertia, w the particle angular speed, F p /M p the force/moment of particle-particle collisions and F w /M w the force/moment of particle-wall collisions.This method enables the dynamic determination of many variables, such as the collision frequencies, the impact velocities, the pressure, the temperature and the RTD signal, but the computation of each particle requires a high computational effort [128][129][130].
In hot-melt extrusion, influences between some modeled material properties of this section can be observed, and recent works study a coupling PBM-DEM to create an accurate modeling [119,121,[131][132][133][134][135]. Figure 23 presents an example (see [136]) of a good combination between these two methods to determine particle physical characteristics.Moreover, computational fluid dynamics (CFD) can also be coupled to PBM-DEM to model particle-fluid interactions in fluidized bed granulation [137][138][139].However, conversely to the previously-introduced methods, PBM-DEM describes only the variations of particle physical properties (assuming the whole device in the steady state) as the granule size distribution (GSD), the moisture content, etc. [119], and demands huge computational efforts, which are not suitable for parameter estimation and control [124].Therefore, related simulation studies are not considered in this work.

Data-Driven Models
In previous sections, most of the models rely on first principles for their derivation.Depending on their degree of details, they range from "grey-box" to "white-box" models.Conversely, "black box" modeling aims at representing the process by its input-output behavior, with no consideration of the system internal state (see Figure 24).In the following, two distinct model classes are presented: • Transfer function models: The input-output behavior of the process around a specific operating point can be described by a transfer function of the form: where Y(s) and U(s) are the Laplace transforms of y(t) and u(t), while G(s) is defined by a set of parameters (a n ...a 0 ; b m ...b 0 ) with n ≥ m.
There exist many techniques to build transfer function models based on output responses to input solicitations, and some of them are summarized in [140].Examples of single-screw extruder modeling can also be found in [141][142][143] and serve the development of more specific techniques related to twin-screw configurations, as in [19,20,144,145].
Common features of these works include: parameter identification based on input step changes; -first-or second-order transfer functions; -inclusion of time delays in some input-output relations (feed flow/die pressure, moisture content/die material temperature, etc).
In [146], a comparative study of different parameter identification methods is achieved.Different algorithms are analyzed in the context of cooking extrusion.The couple "screw rotation speed/motor load" supports relay feedback as a well-performing online method, while batch output error parameters provide an accurate modeling across a large range of operating points.
The main advantages of transfer function models are their simplicity and the availability of a broad array of well-established control techniques, including PID control.However, transfer function models assume that the process behavior is linear in a limited range of operation.If several distinct operating points have to be considered, different transfer functions will usually have to be identified for each of them.In this latter case, gain-scheduling or adaptive controllers will be required.

• Neural networks (NN):
Standard feedforward NNs basically define a nonlinear static map between a selected number of inputs (u) and outputs (y): where i represents a discrete element of sample set I.
One of the most common NN architectures in system modeling is the perceptron [147] consisting of an on/off static function (called the activation function or decision function) delivering a binary output (e.g., zero or one).The sum of a weighted input linear combination is compared to a threshold separating the activation and inactivation zones as in: Perceptron networks may be built using multiple layer structures where all of the neuron outputs depend only on the inputs from the previous layer and do not interact with the same-layer neurons.These structures are called multilayer perceptrons (MLP; [148,149]; Figure 25).The nonlinearity used in the corresponding activation function is continuous (for instance, sigmoid or Gaussian functions).A significant advantage of multilayer NN models is their ability to exploit information contained in every available measurement from nonlinear processes, while, however, the main detrimental consequences are the huge required amount of data to reach acceptable training performances (i.e., the ability to reproduce correctly the outputs) and the fast increasing number of parameters following the addition of layers.The first and the last layers, respectively called the input and output (visible) layers, are distinguished from the intermediate ones, also called the hidden layers.Many other neural network structures, static or dynamic (i.e., using temporally-shifted inputs and outputs), exist, such as the radial basis function network, which, for instance, has proved quite useful in modeling bioprocesses [150] or perceptron-based structures applied to pattern recognition [151,152].These structures generally differ from each other by the activation principle and the training rules.
Steady progress in computational technology allows neural networks to be implemented on-line on many experimental plants [153,154].SISO (single input single output) or MIMO (multi inputs multi outputs) NNs can be used to describe input-output relationships [155][156][157].
In [158], a relation between classical input variables (screw speed, feed flow rate, feed moisture, moisture content, die diameter and temperature) and outlet physical and chemical product properties (radial expansion, density, bulk density, water adsorption and solubility indexes) is established.The selected black-box model is a static three-layered feed-forward network with logistic activation functions trained by a backpropagation algorithm [148].Results show that the method is performing well in reproducing specific material properties and could aim at obtaining generalized models predicting required input parameters to get the desired product properties.
In the following, an illustrative example of the use of transfer function models is given.Using the responses to step changes in the screw rotation speed (Figure 20) and in the input flow rate (Figure 21), transfer functions describing the responses of the die pressure P end and the outlet flow rate Q out are identified using the function "Ident" of the MATLAB identification toolbox.Results are shown in Figures 26 and 27   Second-order transfer functions plus delay are usually sufficient to get an acceptable fit of the input-output behavior.However, as mentioned above, the use of transfer functions is limited to a range of operating conditions where the assumption of linearity is satisfied.Moreover, black-box models only offer an input-output map of the system and no insight into the internal state.

Conclusion
Mathematical modeling of extrusion processes has led to essentially five types of models: • explicit solutions of the mass balance equations, i.e., functions of time describing the residence time distribution; • systems of differential equations, resulting from mass and energy balances expressed around a representation of the systems in the form of continuous-stirred reactors in series; • systems of partial differential equations resulting from mass and energy (possibly momentum) balances in a distributed parameter representation of the system; • PBM-DEM methods, including material properties in differential equation systems describing each particle characteristic evolutions using Newton's second law of motion and the angular momentum • black-box representations of the system, e.g., transfer functions or neural networks linking input and output operational variables, inferred from experimental data and using little or no a priori physical knowledge of the system.
The first class is essentially descriptive, and only the models in the second, third and last categories are of interest in terms of estimation and control.The second class of models can be seen as a first-order finite volume discretization of the partial differential equations involved in the third class of models.Changing the number of CSTRs directly influences both model formulation and the resolution of the solution profiles (numerical dispersion is influenced by the number of discretization cells and can be used in some way to represent physical dispersion).Distributed parameter models offer the advantage of uncoupling model formulation and the numerical solution technique.
Any numerical schemes can be used, ranging from simple finite difference to sophisticated finite element or volume methods.
The PDE models can be complemented by ODEs representing the evolution of moving interfaces, and different models can be used for several sections of the process.The usual spatial representation is in the 1D axial direction, but could be extended to 2D or 3D representations, at the price of increased computational expense and less compatibility with model-based process control strategies, which exploit the simplicity of the model structure and/or sufficiently short computation times.
PBM-DEM methods classify particles with respect to their physical properties (such as granule size, moisture content, porosity, etc.).Differential equations obtained while applying PBM generally require the use of empirical or semi-empirical expressions, functions of particle mechanical characteristics (shear, velocities, etc) determined by DEM.The resulting detrimental effect of this particle-per-particle accurate description is a huge computational effort, inappropriate to control design.
Conversely, data-driven models can be very practical, since they are well-adapted to control frameworks.However, they do not exploit a priori knowledge and offer little insight into the process behavior, outside of the mere input-output dynamics.
To conclude, we primarily recommend a modeling approach based on distributed parameter models.The basic PDEs can be found in [38].These equations together with boundary/continuity conditions can be used to represent several sections of an extruder with different screw configurations.Moreover, the kneading block behavior can be described using assumptions formulated in [55].Finally, the PDEs can be completed by ODEs to determine the spatial extension of fully filled zones and their moving boundaries.The approach proposed in [159] can facilitate the computation.A global model scheme is shown in Figure 28.

Figure 1 .
Figure 1.Sketch of a simple extruder: variable definition.

Figure 2 .
Figure 2. A simple case study.

Figure 3 .
Figure 3. Measurement of the residence time distribution (RTD) at N = 100 RPM and Q in = 0.358 kg/h.

Figure 4 .
Figure 4. Plug flow reactor, continuous stirred tank reactor and actual reactor residence time distributions.

Figure 5 .
Figure 5. Modeling of the residence time distribution through cascades of continuous stirred-tank reactors (CSTRs) or cascades of plug-flow reactor and CSTRs.

Figure 6 .
Figure 6.Discretization of the extrusion device in volume elements, subdivided in cascades of ideal CSTRs.

Figure 7 .
Figure 7. Plug flow reactor and cascade of CSTRs with forward and backward flows.

Figure 8 .
Figure 8. Plug flow reactor in series with a cascade of CSTRs with dead zones.

Figure 10 .
Figure 10.Geometry of a C-chamber.
):the flight gap (with flow Q f ) between the barrel and the screw flight; -the tetrahedron gap (Q t ) between the flight walls; -the calender gap (Q n ) between the flight of one screw and the bottom of the other screw channel; -the side gap (Q s ) between the flanks of the two screws flight; -the channel gap (Q c ) down-channel flow in the screw channel.

Figure 11 .
Figure 11.(a) Five leakage locations in a twin screw extruder; (b) Flow network in a completely filled zone along one pitch length for a single-lobed screw (rectangular boxes represent C-shaped elements; square boxes represent inter-meshing elements).

Figure 12 .
Figure 12.Division of the twin screw in several perfectly-mixed reactors.

Figure 14 .
Figure 14.Die pressure P end and outlet flow Q out evolution following a step change in the screw rotation speed.

Figure 15 .
Figure 15.Die pressure P end and outlet flow Q out evolution following a step change in the feed flow.

Figure 16 .
Figure 16.Representation of the extruder zones.

Figure 18 .
Figure 18.Change of spatial variable χ over the two extruder zones.

Figure 20 .
Figure 20.Die pressure P end and outlet flow Q out evolutions following a step change in the screw rotation speed.

Figure 21 .
Figure 21.Die pressure P end and outlet flow Q out evolutions following a step change in the feed flow step.

Figure 22 .
Figure 22.Scheme of the normal and the tangential contact forces.

Figure 26 .
Figure 26.Die pressure P end and outlet flow rate Q out evolution following a step change in the screw rotation speed.Solid line: transfer function model; dashed line: experimental data.

Figure 27 .
Figure 27.Die pressure P end and outlet flow rate Q out evolution following a step change in the input flow rate.Solid line: transfer function model; dashed line: experimental data.

Figure 28 .
Figure 28.Division of the extrusion system in different sections and moving boundaries of the filled zones.

Table 1 .
Descriptive parameters of the case study.

Table 2 .
Parameter values of the Kumar model.

Table 5 .
Estimated parameter values of the partial differential equation (PDE) model.