Abstract
The present paper deals with the poly(lactic acid) (PLA) reactive processing simulation in a uniform corotating twin-screw extruder that can be readily turned into practical applications in pilot and industrial equipment. The simulator provides a cause-effect guide that can be useful for starting an experimental setup in a reactive screw extruder for a biopolymer in a growing industry. The proposed model considers a free radical ring-opening mechanism involving the main characteristic flows inside the extruder and the non-Newtonian behavior of PLA. The characteristic behavior relating reaction rate, average molecular weights, and polydispersity against chamber number are described by S-shaped and monotonically decreasing curves, for the equipment. Numerical predictions show that this simple and easy to implement model accurately reproduces previously reported data and that the impurity concentration exhibits a marked effect over all the variables, except conversion.
1. Introduction
Poly(lactic acid) (PLA) belongs to the family of aliphatic polyesters, commonly made from α-hydroxy acids, which include polyglycolic or polymandelic acids, and is considered biodegradable and compostable [1]. PLA is a thermoplastic, high-strength, high-modulus polymer that can be made from renewable resources yielding articles for use in either the industrial packaging field or the biocompatible/bioabsorbable medical device market [2]. For further applications, the property profile and price of PLA can be changed by combining PLA with other biocompatible or bioacceptable polymers, fillers, or reinforcements [3, 4]. Alternatively, PLA can be modified by adding plasticisers to obtain more flexible materials [5, 6].
To process PLA on large-scale production lines in applications such as injection molding, blow molding, thermoforming, and extrusion, the polymer must possess adequate thermal stability to prevent degradation and maintain its molecular weight and properties [7]. PLA can be prepared in two ways: either by a polycondensation or by free radical polymerization [7, 8]. The free radicals pathway offers high molecular weight polymers with nearly no remaining monomer. The catalyst currently used industrially is stannous octoate [8], which leads to the ring-opening polymerization of L-lactide and obtains a linear high molecular weight polymer [9]. The choice of initiator system, coinitiator as chain control agent, catalyst concentration, monomer-to-initiator ratio, and polymerization temperature and time, significantly affects the polymer properties, such as the molecular weight, degree of crystallinity, and residual monomer content, which in turn affect the physical-mechanical properties and range of temperature use of the polylactide and its copolymers [10, 11]. Mehta et al. [8] provided a selective summary of the literature on PLA synthesis. The determination of individual rate constants requires the propagation rate constant knowledge, a difficult task that has not been performed properly for this polymer [8]. Evaluation of the rate constants by modeling and simulation in conjunction with the experimental results offers several advantages: initiation, propagation, and termination rate constants can be evaluated/corroborated; the termination mechanism can be ascertained and lastly chain dependent rate constants can be estimated.
One widely used model to predict the complex phenomena inside a twin-screw extruder is the one by Vergnes et al. [12]. They describe the flow for any given screw profile, into a twin-screw extruder, indicating important applications, such as screw profile design, scale-up, compounding, or reactive extrusions.
Reactive extrusion (REX), as the term suggests, involves the synthesis of materials by a melt-phase reaction in an extruder. A detailed summary can be found elsewhere [13–16]. Although reactive extrusion has been established for many years, the use of this technology is currently an active research area [16]. To make PLA production economically viable, the use of a continuous single-step process is highly desirable. The complexity of the flow in twin-screw extruders, as well as the large number of parameters and interrelated variables that affect the flow, makes this process difficult to understand and control. Amongst the underlying flow mechanisms, mixing and residence time distribution (RTD) functions remain preeminent to understand mass transfer. Two broad approaches have been commonly applied to describe the fluid flow mixing. One calculates the velocity and pressure profiles along the screw [12, 16–21], by solving the momentum balance equations. Another approach describes the material flow and mixing by classical flow models widely used in chemical engineering. The parameters of this approach are usually estimated from the residence time distribution (RTD) data obtained from tracer experiments, based on ideal reactor and/or statistical theories [22–24]. Numerical simulation and experimental validations of the local RTD in mixing sections are still rare due to experimental difficulties [22]. Most twin-screw extruders are modular; they are composed of small parts of a given length. These elements are primarily of two types: (i) screw elements with right- or left-handed flights and (ii) kneading elements. In the case of a right-handed element, the screw drags the flow in the downstream direction. The left-handed elements and the kneading elements are restrictive, which tend to hinder the flow of the polymer towards the die, produce complex flows during the rotation, and allow efficient mixing [12, 17]. In REX there are many studies considering isothermal and nonisothermal flow. For the nonisothermal flow the main difficulty arises from the complexity of the velocity field and the streamlines, and thus from the convective term estimation in the heat balance. In fact, the only way to accurately handle this problem consists of considering either a full 3D temperature field or a simplified 1D mean temperature and correcting the temperature dependent variable during the process [12–14, 17]. To simplify the calculations, the isothermal approach is followed in this work.
Banu et al. [25] developed a mathematical model for the L-lactide polymerization process in a twin-screw extruder to investigate both the kinetics and the synthesis of PLA by reactive extrusion. They estimated the RTD during the polymerization to describe the flow using a software called Ludovic [12]. To reduce the cost of manufacturing of PLA, Jacobsen et al. [7] developed a continuous one-stage process using a corotating extruder. This technique requires that the bulk polymerization must be close to completeness within a very short time (5~7 min), which is predetermined by the residence time of the extrusion system.
In this work, the PLA processing simulation in a corotating twin-screw extruder, using a well-known model [12], is performed coupled with a free radical ring-opening kinetic scheme that considers impurities. Numerical predictions are compared and discussed in detail with available experimental data and polymer properties, in terms of the main input variables. The goal of this work is to present a simple-to-implement model for costs and time saving and to analyze further alternatives to improve this attractive sustainable process.
2. Reaction Kinetics
Considerable work has been done on the study of reaction kinetics [6–8, 26]. The nature of the chain-growth process in ring-opening polymerization bears a superficial resemblance to chain polymerization. In the latter, only monomer adds to the growing chain in the propagation step and species larger than monomer do not react with the growing chains. The ring-opening kinetic scheme can be represented as follows.
Initiation is:
Propagation is:
Termination by transfer to monomer is:In (1)–(3), is the initiator concentration, denotes the concentration of polymers of chain length , and denotes monomer concentration. Here , , and are the initiation, propagation, and termination by transfer to monomer kinetic rate constants, respectively. In (3), it is assumed that the charged rings spontaneously form , and is the deactivated polymer with repeat units, which does not participate in any reaction.
Mass balance equations for a batch reactor for this kinetic scheme are: where and are the initial molar concentrations of the initiator and impurities, respectively. According to Kowalski et al. [27, 28] the main impurities are water and octanoic acid. Commercially, SnOct2 contains about 4.5% wt of octanoic acid and up to 0.5 wt% of water. Applying the moment technique to the mass balance equations (4)–(9), and assuming that the rate constants are chain length independent [25], the moment equations are obtained to further evaluate the molecular weights for PLA:In the above equations, (or ) are the live (or dead) th moments and refers to the impurities, respectively. Number () and weight average () molecular weights are calculated as ratios of moments of the molecular weight distribution as follows:where is the monomer molecular weight. Equations (10)–(18) are solved numerically using an iterative procedure described elsewhere [25].
3. Reactive Extrusion
3.1. Flow Assumptions
In general, this model for a corotating twin-screw extruder is well known and has been employed widely for different reaction systems [12, 17, 18]. The main assumptions are (i) isothermal flow, (ii) perfect mixing, (iii) the geometry of the screw which is constant from feed end to tip, and (iv) flow through C-shaped chambers that are modeled as rectangular channels [12].
Partially filled chambers, residence time, and viscous dissipation assumptions are discussed in Sections 3.2.3 and 3.2.4.
3.2. Corotating Twin-Screw Model
In this work, a simplified reactive extrusion model [17, 29] is used. The interaction between the screws is shown in Figure 1.
(a)
(b)
3.2.1. Velocities and Flow Rate Equations
It is well known that an 8-shaped pattern follows the flow along the screws and can be visualized as a succession of elements involving flow through C-shaped chambers and flow through an intermeshing zone between the screws. The flow equations are derived using a cylindrical coordinate system. It is assumed that the flow is locally Newtonian, steady, and isothermal. The radial components of velocity and inertia are assumed negligible. Therefore, the angular (down-channel velocity) and cross-channel velocity components are dependent only on position in the radial direction. Based on these assumptions and their corresponding boundary conditions, the Stokes equations for the down-channel velocity and cross-channel velocity result in the following expressions [12, 17]:where is the viscosity, and correspond to the internal and external radii, respectively, is the radius in a specific point, and and are the axial and tangential pressure drops along the extruder, respectively. Also, and are the relative barrel velocity in the angular and axial directions, respectively. The flow rate through a C-shaped chamber defined in (23) is obtained after integrating (21) with respect to the channel length (). The cross-channel flow rate is equal to zero, the relative velocity of the barrel has the same direction as the main flow, and the pressure change is generally positive:where is the screw rotation speed, is the screw pitch angle, is the angle of a chamber, and and are drag and pressure shape factors, respectively. Tayeb et al. [19, 20] studied the geometry and kinematics in the intermeshing zone using a simple model in which the relative displacement of the screws is neglected and a pressure flow is assumed. Their model uses an approximate mean cross-section defined by the main characteristic parameters of the screw and the barrel. A detailed analysis of the geometry of twin-screw extruders can be found in [12, 18, 19]. The width of this area is calculated using:where (m) is the channel width, is the intermeshed area width, and is the flight angle. The shear rate is calculated using (25), where and are given by (21) and (22), respectively. In these equations, the Newtonian viscosity is replaced by the shear rate dependent viscosity. For this work, the Carreau-Yasuda model is used. Here, the radial velocity is neglected, and the angular and axial velocities are assumed to depend only on the radial position as:It is possible to consider kneading blocks in the present model by using a simplified approach proposed by Vergnes et al. [12], where an equivalent geometry is defined for kneading blocks to simplify the calculations.
3.2.2. Species Balance
The component balance around the th reactor (see Figure 1), assuming equal inlet and outlet flows, is given by:where represents the species in the th chamber and is the polymerization rate.
3.2.3. Residence Time
The chambers in the screw of an extruder can be either partially filled or fully filled. For fully filled subelements, the mean residence time, , is defined as [12, 17]:where is the available volume of the chamber. In partially filled chambers, the reaction mixture is conveyed through the extruder by pure drag flow. Therefore, the mean residence time for an axial length is given by:where is the axial length of the subelement and is the mean velocity along the channel direction. The filling ratio defined as the ratio of the volume occupied by the polymer blend to the free volume available in the chamber is: where is the rotational speed of the screw (given in terms of the number of turns per second) and is the screw pitch.
3.2.4. Density and Viscosity Evaluation
Fang and Hanna [30] studied the rheological properties of amorphous and semicrystalline poly(lactic acid) using a capillary rheometer on an extruder [31]. They analyzed two types of PLA resins (amorphous and semicrystalline) at 150°C and 170°C and different shear rates (30, 50, 70, 90, 110, 130, and 150 rpm screw speed). The viscosity was calculated from the pressure profiles and the volumetric flow rate. They also investigated the melt viscosity as a function of resin type, temperature, and shear rate. Under the same processing conditions, semicrystalline PLA had a higher shear viscosity than the amorphous one. Also, as temperature increased, the shear viscosity decreased for both types of PLA. The PLA melt was characterized as a pseudoplastic, non-Newtonian fluid. Power-law equations were used to describe the behavior of the PLA melts. The authors presented two equations that related the density of the PLA melt to the volumetric flow rate and the PLA melt temperature, . These equations are as follows.
The density is:Banu et al. [25] fitted the experimental data of Cooper-White and Mackay [32] on Carreau-Yasuda’s viscosity model, which is used in this work:The parameters of the viscosity equation (, , and and ) are given in Table 1. Figure 2 shows a plot of the viscosity as a function of the screw speed, where the shear thinning behavior is observed.
3.2.5. Numerical Solution
An iterative procedure is used to solve the set of equations and the calculation sequence is shown in Figure 3. Once the geometric parameters and process conditions are defined, the value of the density is calculated, according to (30), the viscosity at zero shear-rate is considered as initial value (), the monomer, initiator, and impurities concentrations are provided, and all the moment equations (12)–(17) are set to zero. Afterwards, the velocities and flow rate are calculated with (21), (22), and (23), respectively. According to the process conditions, the shear-rate, the filling ratio, and the residence time are evaluated using (25), (29), and (28), respectively. The kinetic equations are solved using (26), and finally the number () and weight average () molecular weights are calculated with (19) and (20), respectively. Figure 3 shows a schematic diagram with the iteration procedure, which is repeated until the value of two successive runs satisfies a tolerance of 1 × 10−7 on the dependent variables.
As depicted in Figure 3, this approach can be applied to other systems due to the fact that any module can be interchanged as required by different process conditions, geometry, momentum equations, rheological model, kinetic equations, residence time, and so on. Therefore, the global model is highly versatile; any module can be modified or completely removed according to a specific process.
4. Results and Discussion
In this section, the numerical predictions for the PLA polymerization reaction in a closely intermeshed corotating twin-screw extruder are presented. The extruder has uniform geometry from feed end to tip with 180 CSTRs or chambers.
4.1. Comparison of the Corotating Screw Extruder Model with Experimental and Numerical Data
Jacobsen et al. [7] performed experiments for the ring-opening polymerization of L-lactide using an equimolar complex of 2-ethylhexanoic acid tin(II) salt, Sn(Oct)2, and triphenylphosphine as catalyst on a corotating, closely intermeshed twin-screw extruder. Their experiments included a specially designed screw concept and catalyst that provided a polymerization propagation rate that was fast enough to use reactive extrusion technology. These authors obtained a number average molecular weight () of 91,100 g/mol and a polydispersity index within a very short residence time (~7 min). However, they did not provide information regarding the value of the kinetic constants. Banu et al. [25] developed a mathematical model for the L-lactide polymerization process in a corotating twin-screw extruder and carried out experimental studies to investigate the process kinetics and the performance of the PLA synthesis process by reactive extrusion. Berzin and Vergnes [29] estimated the RTD of the L-lactide polymerization system which describes the flow in a modular corotating twin-screw extruder. A mathematical model for the isothermal reactive extrusion process was further developed based on the kinetics and RTD, describing the evolutions of monomer conversion and average molecular weight along the extruder, summarized in Tables 2 and 3.
The proposed model is validated by comparing the experimental to predicted monomer conversion and available in the literature [7, 25], using the same kinetic scheme, mass flow rate, screw speed, and monomer/initiator ratio, summarized in Table 2. The values of the rate constant are presented in Table 4.
The obtained values for the monomer conversion, polydispersity index, and and , along the chambers of the extruder, are shown in Figure 4 for two different values for the ratio / (2250 and 4500) using a mass flow rate of 0.75 kg/h and a screw speed of 50 rpm. Clearly, as the ratio / increases, , , and the polydispersity index also increase, while an opposite trend is observed for the monomer conversion. This behavior is the same as the experimental and predicted information presented in Tables 2 and 3. Table 2 compares the predicted conversion obtained by (i) the axial dispersion model; (ii) Ludovic, a commercial software; and (iii) the present model, with the experimental values at two Mo/Io ratios. In the case of this work, where uniform screw geometry is considered, the predicted conversion is very close to the experimental value. It is important to mention that the axial model is simple but requires the information of RTD, and this RTD was obtained from Ludovic, which has a complex geometry screw. Then the axial model is dependent on precise RTD for getting good results. The predictions of the three methods are close to the experimental results obtained using a ratio / = 2250. However, a significant difference is observed among the predicted values of the three methods when compared against the experimental results obtained for the ratio / = 4500. A similar trend is observed for the weight average molecular number when the ratio / = 2250 is increased from 2250 to 4500. This discrepancy is ascribed to the fact that the residence time of the corotating twin-screw extruder configuration used by Banu et al. [25] is larger than that for the constant screw geometry considered in this work (see Table 3). No further attempt was made to fit the data with the present model. However, as explained below, a better fit regarding the molecular weights can be obtained by modifying the impurities content. Another important fact is that viscous dissipation is neglected. This term produces a temperature increase, leading to higher reaction rates.
(a)
(b)
(c)
(d)
4.2. Predictions
4.2.1. Impurities
The modeling of polymerization reactors faces difficulties regarding the description of the flow and mixing, correlation of end-use and molecular properties of the polymer product, and accounting for the role of impurities that often exert a strong influence on the process [25]. Considerable experimental work has been done by various investigators who studied the reaction kinetics of ring-opening polymerization. In experimental studies, various catalyst systems, solvents, and reaction temperature have been used to obtain polymers of molecular weights ranging from a few thousands to over a million. The most frequently used initiator for the L-lactide polymerization process is stannous octoate, its advantages being the high catalytic activity, racemization-free polymerization of L-lactide, and the formation of high molecular weight polymers with high yields [7, 8, 25]. Polymerization with an initial concentration of stannous octoate of 0.05 mol/L with no alcohol added was proved to be very slow by Kowalski et al. [27, 28] because it is certainly initiated by compounds containing hydroxyl groups that are present in the system as impurities. The commercially available stannous octoate contained about 4.5% wt of octanoic acid and up to 0.5 wt% of water. According to the kinetic scheme ((1)–(9)), the termination reactions occur by chain transfer to monomer and to the impurities present in the reaction mixture. However, the impurities do not initiate a new chain acting as an inhibitor. Figure 5 exhibits the impurity concentration effect on the main variables along the extruder. Notice that the conversion remains unaffected in the presence of impurities. As can be seen in (4), the monomer disappearance does not involve the impurities concentration (). Therefore, no effect will be seen in this variable with the increase of . On the other hand, in (8), the dead polymer () increases with the amount of impurities, stopping polymer growth, hence diminishing the molecular weight. Thus, both average molecular weights have an important reduction as the impurity concentration increases, also producing a polymer with a higher polydispersity index.
(a)
(b)
(c)
(d)
4.2.2. Mass Flow Rate
The influence of the mass flow rate through the extruders is analyzed using constant parameter values, namely, = 2250, = 190°C, and the screw speed was set to 50 rpm. Figure 6 shows a reduction in conversion, in , , and the polydispersity index as the mass flow increases. Jacobsen et al. [7] explained this behavior by the effective residence time of the material inside the extruder. A further increase in the throughput would result in a dramatic decrease of conversion, as the viscosity of the resulting material would be too low to provide a successful pressure buildup, which in turn would result in an even lower conversion, leading to a complete breakdown of the reaction. All these experimental evidences are involved and explained semiquantitatively in these reactive extrusion models. According to the reported experiments [7, 25] where the mass flow is increased, first there is a reduction and then an increase in the polydispersity index. They explained that the influence of the intermolecular transesterification reactions, the catalyst, and the characterization technique can generate different behavior in the polydispersity index. The above results show that this simple and easy to implement model, as compared with an expensive commercial software, can be a cost-and-time saving device, helpful to analyze further alternatives to improve this attractive sustainable process. This can be acknowledged by comparing the simulations with the scarce experimental data available [25].
(a)
(b)
(c)
(d)
4.2.3. Screw Speed Impact
According to the experimental results by Jacobsen et al. [7], there is a small reduction in conversion and in the weight average molecular number, as the screw velocity increases (rpm = 50, 100, and 200). A similar behavior was reported by Banu et al. [25]. This behavior is predicted also in this work; however, the screw speed has a more drastic effect in conversion, as shown in Figure 7. When increasing the rpms, the volumetric flow augments also, and the degree of conversion is reduced. This fact is explained by the effective residence time of the material inside the extruder [7]. Equations (27) and (28) show how the residence time decreases with increasing throughput and increasing screw speed, respectively. If the chambers are modeled as partially filled (equation (28)), then the possibility to have more than one entry along the process exists. Then, the filling ratio in the starved elements can be calculated using (equation (28)). Decreasing the residence time yields an incomplete polymerization and also lower and . With an increment of the screw speed, a low reduction on occurred experimentally. According to Jacobsen et al. [7] the reduction in is due to the heat generated by viscous dissipation, which is also a source of energy for further enhancement of the propagation step. Our model prediction underestimates the average molecular weight. However, the qualitative trend is correct. The polydispersity index is also reduced when rpms increase, as shown in Figure 7.
(a)
(b)
(c)
(d)
5. Conclusions
The simulation for the PLA production in a corotating twin-screw extruder has been performed. The corotating extruder model is compared with experimental and numerical results from the literature, and a good match is found with the scarce existing experimental evidence. The influences of the main variables, such as the screw speed, mass flow rate, / ratio, and impurities, have been assessed. The numerical results verify that the PLA polymerization can be performed under a free radical ring-opening polymerization scheme for a reactive extrusion process. According to the experimental data reported in the literature, the catalyst is very important for the reaction time, allowing high conversion and high molecular weights. Results show that the screw speed is an important parameter for reactive extrusion process. High rpms allow working at lower power consumption, but this decreases the residence time and lowers conversion and thus . The use of the Carreau-Yasuda viscosity model for viscosity evaluation was considered because it involves the evolution of the reactive process. Increasing the flow rate produces a reduction in conversion and , due to a reduction in the residence time. The commercially available stannous octoate contains low concentrations of impurities which reduce the molecular weight drastically. These cause-effect relationships presented here should be useful as a guide for starting an experimental setup in a reactive screw extruder. There are many catalysts which have been reported to initiate this polymerization. Many current PLA polymerization methods employ stannous octoate in different solvents as the catalyst, and it has been shown to be very effective. However, there is still a need for a safer, faster, and more stable catalyst and the evaluation of the individual kinetic rate constants that allow for better model predictions. The comparison with the scarce experimental data shows that this simple and easy to implement model can be a cost-and-time saving device. The extruder configuration captures the physics of the problem and, therefore, the operating conditions can be consistently manipulated to reduce experimental runs, avoiding waste materials and reducing operational costs.
List of Symbols
| Parameter of the Carreau-Yasuda model | |
| : | Pitch of the screw |
| CSTR: | Continuous stirred tank reactor |
| : | Heat capacity, J Kg−1K−1 |
| : | Screw diameter, m |
| : | Energy dissipation, J s−1 |
| : | Shape factor for the drag term |
| : | Shape factor for the pressure term |
| : | Height of the chamber, m |
| : | Initiator concentration, mol L−1 |
| : | Initial initiator concentration, mol L−1 |
| : | Kinetic rate constant for initiator, L mol−1 min−1 |
| : | Kinetic rate constant for propagation, L mol−1 min−1 |
| : | Kinetic rate constant for chain transfer to monomer, L mol−1 min−1 |
| : | Kinetic rate constant for chain transfer to impurities, L mol−1 min−1 |
| : | Length in the -direction, m |
| : | Mass flow rate, kg h−1 |
| : | Parameter of the Carreau-Yasuda model |
| : | Monomer concentration, mol L−1 |
| : | Initial monomer concentration, mol L−1 |
| : | Concentration of deactivated polymer with repeat units, mol L−1 |
| : | Monomer molecular weight, kg mol−1 |
| : | Number-average molecular weight, kg mol−1 |
| : | Weight-average molecular weight, kg mol−1 |
| : | Power-law constant |
| : | Rotational rate, min−1 |
| PLA: | Poly(lactic acid) |
| : | Concentration of activated monomers/radicals, mol L−1 |
| : | Concentration of polymers of chain length , mol L−1 |
| : | Tangential pressure drop |
| : | Axial pressure drop |
| : | Volumetric flow rate, m3 h−1 |
| : | Flow through a C-shaped chamber, m3/h |
| REX: | Reactive extrusion |
| RTD: | Residence the time distribution |
| : | Radius, m |
| : | Screw external radius, m |
| : | Screw internal radius, m |
| : | Polymerization rate, mol L−1 s−1 |
| : | Impurities concentration, mol L−1 |
| : | Initial impurities concentration, mol L−1 |
| : | Temperature of the melt, °C |
| : | Adiabatic temperature rise, °C |
| : | Mean residence time, h |
| : | Mean velocity along the channel direction, m h−1 |
| : | Down-channel velocity, m h−1 |
| : | Chamber volume, m3 |
| : | Relative barrel velocity in the angular direction, m h−1 |
| : | Cross-channel velocity, m h−1 |
| : | Width of the channel, m |
| : | Width of the intermeshed area, m |
| : | Relative barrel velocity in the axial direction, m h−1 |
| : | Required power per unit volume to rotate the screw |
| : | Mass of species in the th chamber |
| : | Channel length, m. |
| : | Flight angle (degrees) |
| : | Index |
| : | Shear rate, s−1 |
| : | Parameter of the Carreau-Yasuda model |
| : | th moment of live polymer radical distribution, mol L−1 |
| : | th moments of dead polymer radical distribution, mol L−1 |
| : | Non-Newtonian viscosity, Pa s |
| : | Zero-shear viscosity, Pa s |
| : | Screw pitch angle |
| : | Density, kg m−3 |
| : | Screw rotation speed, rad s−1 |
| : | Angle of a chamber, degrees. |
| : | Initial condition |
| : | th chamber or th CSTR. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work has been supported by the Consejo Nacional de Ciencia y Tecnologia (CONACyT), Grant no. 167474. This project was supported by the research Grant 20152030 SIP-IPN, Mexico. Francisco López-Serrano gratefully acknowledges support by PAPIIT IN113215, UNAM. The authors would like to express their gratitude to Ing. Ruben Mil for his contributions to this work.