DETERMINING FINAL PRODUCT PROPERTIES AND KINETICS STUDIES OF POLYPROPYLENE POLYMERIZATION BY A VALIDATED MATHEMATICAL MODEL

In spite of a long history of using Ziegler-Natta catalyst for producing polypropylene, because of the complexity of the chemical and physical variables and the non-linearity in their effects on the polymerization, nobody has prepared a validated model which able to predict the final product properties and the critical kinetic indices at the same time. In this work, a mathematical model for slurry polymerization was implemented by polymer moment balance technique (population balance approach) in MATLAB/SIMULINK environment. The modeling approach is based on the expected value in statistics. Its advantage is that some important properties such as number average and weight average molecular weights could be directly calculated from kinetics equations. Other significant polymer indices, such as polydispersity and melt flow index which are a function of molecular weight, may be easily calculated and predicted by this model. In addition; by using the reaction profile rates (Rp-t) curves which the model is able to calculate and plot them, it is possible to determine the kinetics indices such as the initial reaction rate, deactivation constant, activation energy, and yield of the catalyst. The model was validated by experimental data obtained from a lab scale semi-batch reactor.


INTRODUCTION
Nowadays, polypropylene is one of the widely used polyolefin products, and its application is strongly affected by the final product properties.In spite of sixty years of history of the use of Ziegler-Natta catalyst in producing polypropylene, the polymerization system performance still remains as a black box [1] due to the fact that the kinetics of polymerization with Ziegler-Natta catalysts are quite complex [2].Like other polyolefins, the most vital indices of polypropylene properties are comprised of melt flow index (MFI); number and weight average molecular weights (Mw and Mn), and polydispersity index (PDI).Since more kinetic data are often obtained from polymerization rate profiles, it is necessary to predict the rate of polymerization by the model.The vital kinetic parameters are the reaction rate-time (Rp-t) curves, initial reaction rate (Rp0), deactivation constant (K d ), activation energy (E a ), and yield (Y).
So far, the majority of researchers, to understand the polymerization system, carried out their research in a conventional manner based on trial and error or experimental approach; though this methodology is not reliable and applicable because the results are heavily dependent on the test and laboratory conditions.Hence, a validated mathematical model that could be capable of predicting the final product properties and the vital kinetic parameters is desirable.
There have been many studies on different aspects of polypropylene polymerization; but so far, only a very few studies have been carried out regarding polyolefin reactor modeling; and most of them focused on heat and mass transfer inside the catalyst particles or the mechanism of reaction, and or how to grow the polymer particles via multigrain; On the other hand, their models have not been validated with experimental data [3][4][5][6][7][8].Some other scholars have solely centered on a loop or fluidized-bed reactors (FBRS), namely bulk or gas phase polymerization [2,[9][10].Since the researchers' mesoscale viewpoint and concern on heat and mass transfer limitations to the catalyst particles, their models were not applicable for the macro-scale model to determine final product properties and also the overall kinetics study of the polymerization.
Because using the profile of polymerization rate, most of the vital kinetics parameters could be obtained; no validated mathematical model is proposed to plot the rate of propylene polymerization which is the aims of this study.Up to now, the profile rate of the polymerization is determined by the experimental study.Samson et al., the first researcher, have shown the rate of the polymerization in the liquid phase experimentally [11].
On the other hand, the effect of hydrogen on the rate of the polymerization and final product properties is significant to be considered; because hydrogen as transfer agent has a direct effect on final properties, the activated site of the catalyst used and also the rate of the polymerization.
Regarding the effect of hydrogen on the kinetics, some experimental investigations have been carried out [12][13][14]; but the results of them have been ambiguous and even contradictory.Guastalla and Gianinni concluded that the initial rate and activity of the catalyst increase about 2.5 times when hydrogen exists in the polymerization reactor [15].Spitz et al. reported that low hydrogen concentrations in the reactor caused to the enhancement of the rate profile, higher hydrogen levels lowered activity and increased deactivation of the catalyst used [16].Rishina et al. confirmed the previous work; however, they pointed out that the role of hydrogen is temporary [17].Contrary, Soga and Siona concluded that propylene polymerization rate decreases with increasing hydrogen partial pressure [18].Kahrman et al. obtained different conclusions that hydrogen had not only an effect upon the polymerization rate for low hydrogen concentrations but also the rate of polymerization decreases at high hydrogen concentration [19].
Al-haj Ali et al. proposed a generalized model for hydrogen response based on the dormant site theory in liquid propylene polymerization.However, his research work was only based on experimental data without providing a mathematical model which might be able to predict the polymerization rate profile and the vital indices of final product properties [14].Reginato et al. modeled an industrial-scale loop reactor by using a nonideal continuous stirred tank model to explain the industrial process and compared their simulation results with commercial plant data [2].The research targets have been defined to predict the macroscopic of the process, dynamics of the plant, advanced control strategies, grade transitions, and average polymer properties.In their paper, despite offering formulas for some important final properties such as number average molecular weight (Mn); weight averages molecular weight (Mw) and polydispersity index (PDI); they did not mention about validating the model.
Varshouee et al. proposed a validated model which was able to determine the effect of the amount of hydrogen on the percent of activated sites on the catalyst used [20].Despite the rather long history of the usage of Ziegler-Natta catalyst for producing polypropylene, the overall performance of the polymerization systems still remains as a black box [1].Hence, unlocking the black box by means of a validated model is the best way.Up to now; no one, to the best of our knowledge, has studied and provided a validated model capable to predict the vital indices of the final product properties and crucial kinetic parameters at the same time.In view of the significance of the matters, with this paper, we solve the prevailing issues via providing a validated model.The selected modeling technique is the polymer moment balance method (population balance approach) coded in MATLAB/SIMULINK for slurry polymerization.By this model, the most vital indices of polypropylene properties such as melt flow index (MFI); number and weight average molecular weights (Mw and Mn), and polydispersity index (PDI) could be calculated directly.In addition, the model computes and plots the reaction rate-time (Rp-t) curves for the Kinetics studies indices, such as the initial reaction rate, deactivation constant, activation energy, and yield of the catalyst.The model was validated by experimental data obtained from a lab scale semi-batch reactor by using the 4th generation of Ziegler-Natta catalyst with an acceptable margin of error.The model was validated through experimental data obtained from a lab scale semi-batch reactor by the usage of the 4th generation of Ziegler-Natta catalyst with an acceptable margin of error.By means of the model, it is easily possible to adjust process conditions for the favored product without hazard.In conclusion, optimum process conditions for maximum yield of catalyst and favored polymer properties can be determined by means of the model.For this study, the optimum condition is 70 o C and 18 mg of hydrogen in the system.

Materials
The materials used in this study are as follows: the catalyst used; the 4th generation of spherical MgCl 2 supported Ziegler-Natta catalyst containing 3.6 wt% Ti.Diisobutyl phthalate (DIBP) as internal donor supplied by Sudchemie, Germany.Triethylaluminium (TEAL of 98% purity) as a cocatalyst from Merck, Germany which diluted in n-heptane.Cyclohexyl methyl dimethoxy silane (CMDS) as an external donor; from Merck, Germany.Polymer-grade propylene was provided from Shazand Petrochemical, Iran.Hydrogen and nitrogen with, >99.999% purity.

Polymer synthesis
In this study, the homopolymerization was carried out in n-heptane media as a slurry.Polymerization reactor was a 1 L stainless steel vessel manufactured by Buchi Co. as shown in Figure 1.A high-pressure N 2 line was used to transfer liquid monomer and catalytic system into the reactor.The catalyst was injected into the reactor through a stainless steel cylinder under N 2 atmosphere.The individual gases were then filtered and flow of each reactant was measured and controlled with a mass flow controller manufactured by Brooks Company.
Rp-t curve data directly come from the monitor of the polymerization set-up.The molecular weight of products are measured by gel permeation chromatography (GPC), Agilent PL-220 model which equipped with TSK columns at 155 o C using 1,2,4-trichlorobenzene as a solvent.MFI of samples is evaluated according to ASTM 1238 at the temperature of 230 ºC and load of 2.16 kg.

Polymerization procedure
A typical polymerization procedure consisted of reactor preparation, polymerization, and discharge steps.Details are as follows: Firstly, the reactor was flushed with nitrogen gas for 1 hour at 90 o C and was reduced reactor temperature until 20 o C, then purged with propylene gas in 15 min.Afterward, 500 mL heptanes as a solvent were introduced to the reactor; next, all the inputs and outputs of the reactor were closed and were stirred at 200 rpm in 5 min for solvent degassing under a vacuum pump.Subsequently, hydrogen was entered to the reactor (based on recipe condition), then propylene was introduced to reactor according to the installed control program, then the reactor was heated up until reach to equilibrium thermodynamic conditions (T = 70 o C, Pr = 7.5 bar), finally the reactor was ready for injecting catalyst to start polymerization.Injecting catalyst to the reactor was carried out under pressure via an injection system during polymerization time (two hours) at constant temperature and pressure, that is to say, the reactor was executed under isothermal and isobaric reactor condition.Data were collected every five seconds.

Assumptions
The following modeling assumptions are considered: (1) It was supposed that propylene polymerization was carried out in the amorphous phase and amorphous phase concentrations During polypropylene polymerization is at the thermodynamic equilibrium condition that obeys from Sanchez and Lacombe Equation (SLE) [21] for calculating the amount of X = CH/Cm, the hydrogen molar ratio.(2) It was assumed that γ1 = γ2 =.....= γnc.Where γ is equilibrium constant and NC is a number of solvent in slurry phase components [2].(3) The reaction temperature, pressure, and monomer concentration were kept constant during the polymerization process.(4) The resistance of both mass and heat transfer and the diffusion effect of the reactants were ignored.( 5) It was assumed that the propagation constant is independent of the length of the growing polymer chain.(6) Using "dormant sites theory" for activating catalyst by hydrogen concentration [11].

Mathematical formulation
Olefin polymerization kinetics with Ziegler-Natta catalysts might be fairly complicated [2].To date, several reaction steps have been proposed in the open literature [2,7].However, the most comprehensive steps were proposed by Zacca [2].The ODE mass balance equations in the model are as follows: Since the model is a semi-batch process and also to be assumed that the monomer concentrations are constant during the polymerization, the input and output terms are eliminated (Qf and Q0), therefore, the terms of η and ζ are meaningless in this study.Table 1 shows possible reactions with their rate equations in the polymerization reactor [2].
Table 1.The probable reactions and their rate equations in propylene polymerization used in the model [2].

Reaction step
Component Reaction Rate equation Site activation Hydrogen

Chain transfer
Hydrogen The component rate equations used in the model are listed in Table 2 as follows [2].
Table 2.The component rate and moment equations used in the model [2].

Moments equations
Live polymer The time-dependent concentration variations used in the modeling are as follows: Where: k is the site number of the catalyst.
In this study, it is supposed that the catalyst has mono-site, and then k is equal to one.Here, CH, CA, CE, CMi, CB, CS, Ccat, and P0 is the concentration of hydrogen, co-catalyst (aluminum alkyl), electron donor, monomer, poison, site transfer, catalyst and potential site in the polymerization in respectively.The component rate equations and moment equations used in the model are listed in Table 2.The final product properties of polypropylene can be estimated by the moment equations.The basic polymer properties, called as end-use properties, are four items; number average molecular weight (Mn), weight average molecular weight (Mw), melt flow index (MFI) and polydispersity index (PDI).The relationship between the moment and these indices are defined by the following equations: The basic polymer properties, called as the final product properties, are four items; number average molecular weight (Mn), weight average molecular weight (Mw), melt flow index (MFI), and polydispersity index (PDI).The Equations used in the model are as follows [2]: used to adjust kinetic parameters [20].

Modeling algorithm
This study has outlined the algorithm for programming the mathematical model in a MATLAB/SIMULINK environment, as shown in Figure 2a.It is composed of two part; mainprogram (as named "Runsim") and subroutine (function file).For obtaining kinetic constants in the model, it is proposed a new approach as iterative method algorithm by using consistency property of ODE's equation in Figure 2b.In this study, From open literature [2,9], the initial guess of kinetic constants was estimated and applied to the model, next the constants were exactly adjusted and determined in accordance with the catalyst used in the set-up (experimental data ) by the proposed algorithm in Figure 2b.Comparing the polymerization profile rate of the model outputs and the experimental data in Figure 3a and 3b implies that the fairly accurate kinetic constants have been adjusted and applied in the model.

Determination of the kinetic parameters; Rp0, K d , E a and Y
As shown in Figure 3-a, 3-b, a typical polymerization rate profile is comprised of two zones; (I) initial polymerization; the startup zone and (II) the quasi-steady-state zone.Each zone has a significant meaning in the kinetic analysis.Detailed discussions on these issues have a considerable debate in the literature and are not repeated here for the sake of brevity [14].Pater et al. and some other researchers have shown that the rate of polymerization at isothermal conditions can be described as a first-order process in monomer concentration (Eq.9) and the deactivation of the catalyst as a first-order process with the number of active sites (Eq.10) [12][13][14]  After integrating and combining the above equations at isothermal conditions leads to the following equation.
) exp( Equation ( 11) is often used to describe the time-dependent rate of polymerization at isothermal conditions; but at non-isothermal, it is necessary to use the following equation well-known to the temperature-dependent rate of polymerization based on Arrhenius equation.
Where Rp0 is the initial reaction rate, K d is the deactivation constant and t shows time in equation (11); and also E a,d is the activation energy for the lumped deactivation reaction, T indicates the temperature and Rp0 is the initial reaction rate in equation ( 12).These parameters can be calculated based on the profile rate curve of polymerization (t = 0).
Since the rate of polymerization has a dependency on temperature, with Rp0 at different temperatures and by using Arrhenius equation, the activation energy of any type of catalyst is easily predictable (Figure 4 for this study) [20].Two very important points should be considered about activation energy of Z.N catalyst; first, it is independent of the hydrogen concentration [14], and second, this is an intrinsic property of any catalyst and thus can be expected to differ from one type to another.

589
The yield of the polymerization can be calculated by integrating the rate.(14) Amount of Y calc is exactly equal to the area under the profile curve.If this value is multiplied by the amount by weight of the catalyst, the amount of polymer produced will be obtained in each batch.In the experimental part, the yield is measured by weighing the dry product from batch polymerization.But in fact, much more monomers are entered to the reactor, called as the consumed monomer, part of them are reacted and the other part remains as an unconverted monomer in liquid and gas phases.The model is capable to give yield and consumed monomer, directly by calculations.By plotting the natural logarithm of the reaction rate versus polymerization time, a linear fit can be made where the slope of the fit line is K d .

RESULTS AND DISCUSSION
In the beginning, it is necessary to be shown that the performance of this model is acceptable and validated in a proper way.Comparing the profile rates of the polymerization come from the model output and from experimental work might be the most effective manner for investigation the performance of the model.A typical comparison is shown in Figure 2(a) and (b), which illustrate the comparisons of the model and experimental profile rate of the polymerization at 70 ºC in the absence and presence of hydrogen; they show that the model trend has an acceptable error and is in line with the experimental curve.Therefore, it is concluded that the model has an acceptable performance and been validated in a proper manner.Existing errors might be explained by the following reasons: (i) truncation, round-off and method errors in numerical calculation, (ii) personal and measurement equipment errors, (iii) the errors inherent in the equation of state selected and (iv) the assumptions error.
For the reason that activation energy is an inherent property of the reaction and is independent of any process variable, it may be considered as another reason for indicating of the validated model.As a result, the activation energy of the model and experimental works are calculated in accordance with Figure 4 and compared with other works in Table 3. Table 3. Reported activation energies (Ea) in propylene polymerization [20].

Worker
Catalyst system Phase Ea (kJ/mol) Yuan et al. [22] δ-TiCl3.l/3 AlCl3/DEAC Slurry 53.9 Soares et.al. [23]  As it can be seen in Table 3, the activation energy which has been calculated from the model and experimental works are in line with the other researchers' reports.This is another reason that the model is validated.The model output and experimental results at different conditions are summarized in Table 4.
Comparing the amount of the model output and experimental results in Table 4 indicate the accuracy and performance of the model.As can be seen from Table 4, the experimental results and model outputs have an acceptable margin of error in regard to the most vital indices of polypropylene properties such as melt flow index (MFI); number and weight average molecular weights (Mw and Mn), and polydispersity index (PDI) which are defined as the main and the first purpose of this work.Consequently, it is concluded that the model is easily able to calculate and predict the vital indices of final product properties.Investigating indices of kinetics, such as the initial reaction rate (Rp0), deactivation constant (K d ), activation energy (E a ), and yield (Y) are the second purpose of this work.As mentioned earlier, the initial reaction rate (Rp0) and deactivation constant (K d ) are obtained by plotting the natural logarithm of the reaction rate versus polymerization time, a linear fit can be made where the slope of the fit line is deactivation constant (K d ) and the intercept of the line is initial reaction rate (Rp0).For each run, in two cases, i.e. model and experimental, the value of reaction rate (Rp0) and deactivation constant (K d ) are separately calculated and inserted in Table 4.According to Equation (14), the amount of Y calc is exactly equal to the area under the profile curve which been calculated and inserted in Table 4 as well.As indicated in Table 4, it is concluded that the model is easily able to calculate and predict the vital indices of kinetics study as well.
Figure 5(a) and (b) present the effect of the reaction temperature and hydrogen amount on the yield of the catalyst used.Figure 5(a) reveals that the yield of the catalyst at 70 ºC has the maximum amount.It means that the reaction temperature at 70 ºC could be considered as the optimum temperature in this work.The effect of hydrogen changes on the yield of the catalyst used at 70 ºC was plotted in Figure 5(b), it implies that the maximum of the yield was obtained at 70 ºC and 18 mg hydrogen which might be considered as an optimum condition in this study.On the other hand, this conclusion is confirmed by evaluating the effect of reaction temperature and hydrogen amount on the deactivation constant (K d ) as the following figure.
Increasing temperature according to Arrhenius equation Equation ( 13) cause to increase deactivation constant (K d ) which no desirable event, because the age of the catalyst used would be reduced.This event is affirmed in Figure 6(a); the figure shows that the maximum of the deactivation constant (K d ) happens at 75 ºC, that is to say, the catalyst has the minimum life in a reactor which is an undesirable event.While at 70 o C it is more moderated and suitable; this might be considered as another reason which the reaction temperature at 70 o C is the optimum temperature.In another word, after the optimum temperature (i.e.70 o C), increasing K d give rise to decrease the yield of polymerization.This event is clearly shown in Figures 5a and 6a with comparing them.The decrease in polymer yield at temperatures higher > 70 o C is due to catalyst deactivation either by over reduction of the catalyst sites or due to alkylation with the Lewis base [11,24].The rate of reaction is proportional to the monomer concentration absorbed by the amorphous phase of the polymer.The monomer concentration in amorphous polymer decreases with increasing temperature, and consequently, the reaction rate and yield decrease if the reaction temperature increases after the optimum temperature (i.e.70 o C) as revealed in Figure 5(a).
Table 4 also shows the effect of hydrogen on the kinetic parameters of the polymerization.By increasing the hydrogen amount, the peak of the initial reaction (Rp0) rate shifted towards higher values.This means that increasing hydrogen amount leads to increase the yield and the rate of polymerization.As shown in Figure 5(b), the maximum yield is attained at the 18 mg hydrogen amount.However, after that, the yield of polymerization decreases.This behavior exactly occurs for K d (Figure 6-b).Therefore, it is concluded that the 18 mg value is an optimum concentration of hydrogen for the polymerization system.
These findings, i.e., the decrease in the rate at high hydrogen concentration, were observed earlier by other workers for gas phase polymerization and slurry polymerization as well [14,25].Natta et al. attributed this decrease in polymerization rate to the slow reactivation of the catalyst-H bond, which results from a chain transfer reaction, or a side reaction such as partial hydrogenation of aluminum alkyl, since the activation energy is virtually independent of the hydrogen concentration [14].Therefore, the calculation of activation energy is not different in the absence or presence of hydrogen.The calculation is based on the Arrhenius plot of the initial polymerization rates Rp0 at different temperatures (as shown by Figure 4), and our results (experimental and model) are compared with the results of other researchers in Table 3.The effects of hydrogen and temperature on final product properties were also investigated.The results inserted in Table 4, demonstrate that Mw has a proportional relation with the reaction temperature and consequently, according to Equation ( 8), MFI decreases by increasing temperature.It is observed that PDI has a minimum amount at 70 o C, namely, at the condition the dispersity polymer chain is optimum.The constants of the Equation ( 8) were calculated in accordance with Figure 7 [20].The results summarized in Table 4 reveals that the molecular weight of polypropylene decreases with increasing hydrogen amount due to a chain transfer agent of hydrogen during olefin polymerization.

CONCLUSION
The present work is designed to investigate final product properties and kinetics studies by using a mathematical model which validated by experimental data obtained from a lab scale isothermal semi-batch reactor by using the 4th generation of Ziegler-Natta catalyst with an acceptable margin of error.The model predictions were in good agreement with the experimental results and revealed that the optimum temperature and hydrogen concentration values were 70 o C and 18 mg, respectively.At the optimum temperature, the PDI was at a minimum amount, indicating optimum dispersity of polymer chains.In the absence of hydrogen, the increase in temperature caused to reduce of molecular weight and enhancement of MFI.Since the energy of activation is independent of hydrogen and can be a criterion for determining the validity of the model.The model is able to compute and predict the indices of kinetics, such as the initial reaction rate (Rp0), deactivation constant (K d ), activation energy (E a ), and yield (Y) and also could to be used for evaluating the performance of new or unknown catalyst in the presence or absence of hydrogen in the polymerization.We suggest an investigation of the catalyst deactivation model for future studies.

Figure 1 .
Figure 1.Simplified schematic of the reactor system.


Determining final product properties and kinetics studies of polypropylene polymerization Bull.Chem.Soc.Ethiop.2018, 32(3) , and QR are feed volumetric flow rate, reactor-output volumetric flow rate and volumetric recirculation flow rate in respectively.In which VR and Rj are defined as reactor volume and j component reaction rate in equation 1.

( 8 )Figure 2 .
Figure 2. (a) The general algorithm modeling in this work and (b) the iterative methodology used to adjust kinetic parameters [20].

Figure 3 .
Figure 3.Comparison of the experimental and model profile rates at 70 o C, (a) in the absence of hydrogen; (b) in the presence of hydrogen. : Gholam Hossain Varshouee et al.Bull.Chem.Soc.Ethiop.2018, 32(3)

Figure 5 .
Figure 5.The comparison of experimental data and model prediction of yield (gr PP): (a) in the absence of hydrogen with the variation of temperature (Runs 1, 2, 3) and (b) at the optimum temperature (70 ºC ) with the variation of hydrogen amount (Runs 2, 4, 5).

Figure 6 .
Figure 6.The comparison of experimental and model deactivation constant, (a) in the absence of hydrogen and (b) in the presence of hydrogen.

Figure 7 .
Figure 7. Fitting melt flow index (MFI) for polypropylene samples in the absence of hydrogen and different temperature for obtaining the constants of Equation 8 (Runs 1, 2, 3).

Table 4 .
The comparison of model output and experimental results at different conditions.