Modeling and Simulation of an Isothermal Suspension Polymerization Reactor for PMMA Production Using Python

Show more

1. Introduction

Poly (methyl methacrylate) (PMMA) is a transparent thermoplastic used in a wide range of fields including medicine, art, and aesthetics. It is also used as an alternative to glass due to its shatter-resistance property [1] . There has been a widespread need to use PMMA based components in Uganda especially for construction and medicine purposes, however advancing the technology that addresses industrial polymerization processes is often hampered by a lack of fundamental research activities in academic and the few industrial labs available in the country.

Industrially, PMMA is obtained from the free-radical polymerization of methyl methacrylate [2] which may be produced via methanolysis of methylcrylamide sulphate [3] . The polymerization reaction is commonly initiated by thermal radical forming agents often employing organic peroxy compounds and azo-compounds [4] . Suspension polymerization is used for the commercial production of many important polymers other than PMMA such as polyvinyl chloride and polystyrene [5] with similar properties to those produced from bulk polymerization. The limitations of increased viscosity, reaction temperature, and volume contraction in the latter process are always avoided [6] . The process may be described as being operated batch-wise however, in essence it is often in a semi-batch mode since some material enters the reactor after the start of polymerization.

During the process, the monomer is stirred with about twice its volume of water and dispersants forming droplet-like distribution of the monomer phase in water. The mechanism of polymerization corresponds to that of bulk polymerization and as such, the droplets become increasingly viscous during the course of polymerization. Distributors such as gelatin stabilize the suspension and prevent the droplets from adhering to one another during collisions. The system is heterogeneous therefore the polymer is collected as granular beads on completion of polymerization and it is washed to remove adhering distributor and salts then dried [7] .

Various works have been published on the modeling of methyl methacrylate polymerization processes. A major feature of these studies is the observation in the increase of the mass reaction viscosity with monomer conversion resulting into a deviation from the “normal” kinetic. The severe reduction in the mobility of the macroradicals auto accelerates the reaction causing an increase in the gel effect and a decrease in the termination rate constant. If the polymerization is carried out at lower temperatures, the glass effect occurs when the transition state is reached at a certain conversion causing a decrease in the propagation rate constant. This results in the interruption of the reaction before the monomer is completely consumed. The increased viscosity at high monomer conversion leads to a decrease in the initiator efficiency due to the cage effect [8] . In most of the reviewed works, gel and glass effects are modeled according to Chiu et al. [9] and the cage effect is modeled according to Achilias-Kiparissides [10] . Most of the models considered are reported to have worked well for MMA polymerization [11] , however the general model developed by Seth and Gupta [12] and reworked by Ghosh et al. [13] using free volume theory is noted as being most superior [8] . Several simulation packages to assess the performance of the models under consideration have been used such as MATLAB [14] , gPROMS [15] , FORTRAN [16] , which are based on utilization of sequential quadratic programming, ordinary differential equations and fuzzy-neural modeling.

From the above, it can be seen that techniques applied in modeling polymerization processes are reasonably well developed and several commercial simulation packages are available. However a user of such techniques and packages should keep in mind that software cannot be a substitute for engineering judgement and its use without understanding is a dangerous abandonment of professional responsibility [17] .

The aim of this paper is to model and simulate the polymerization of methyl methacrylate to PMMA in an isothermal suspension reactor. The reactor performance is studied by modeling and simulation using Python 3.5 instead of using experimental assessment. This is because the procedure requires low computational demand in performing the reactor studies and besides high-quality data is available for validation and comparison with literature values.

2. Model Formulation

The kinetics of MMA polymerization are described basing on the free radical polymerization mechanism scheme in Table 1 below. The scheme map comprises the initiation step, where the initiator (I) decomposes to form reactive radicals (R), the addition of monomer molecules (M) to the reactive polymer chain formed from the reactive radicals and the termination step where deactivation of polymer radicals occur. The initiator used is Benzoyl Peroxide (BPO).

P_{ }and D represent the live and dead polymer with n and m monomeric units respectively; whereas the kinetic constants for the initiator decomposition, initiation, propagation and termination are K_{d}, K_{i}, K_{p} and K_{t} respectively.

In this work, the process model in [12] [13] with slight adjustments was followed. The assumptions used in this work are:

1) The reactor operates isothermally at 70˚C

2) All reaction steps in the system are irreversible

3) Chain transfer to the monomer is negligible compared to other reaction steps

4) The dispersion generated in water isolates the monomer thus negligible evaporation rate

5) Each droplet behaves as a batch reactor operating in bulk

6) Termination is due to disproportionation only

7) The volume of the liquid phase is unaffected by the small amount of initiator

8) There is no polymer at the beginning of the process

Using the above assumptions, the initiator (I), monomer (M) and radical species (R) balances are;

Table 1. Kinetic scheme for the free-radical polymerization of MMA.

$\frac{\text{d}I}{\text{d}t}=-{K}_{d}I+\text{input}$ (1)

$\frac{\text{d}M}{\text{d}t}=-{K}_{p}\left(\frac{{\lambda}_{0}M}{V}\right)-{K}_{i}\left(\frac{RM}{V}\right)$ (2)

$\frac{\text{d}R}{\text{d}t}=2f{K}_{d}I-{K}_{i}\left(\frac{RM}{V}\right)$ (3)

The method of moments is invoked in order to reduce the infinite system of molar balance equations required to describe the molecular weight distribution [18] that would have presented roughly 10,000 - 50,000 stiff differential equations to be solved simultaneously [19] . Therefore the first three radical species and dead polymer moment balances are;

$\frac{\text{d}{\lambda}_{0}}{\text{d}t}={K}_{i}\left(\frac{RM}{V}\right)-{K}_{t}\left(\frac{{\lambda}_{0}^{2}}{V}\right)$ (4)

$\frac{\text{d}{\lambda}_{1}}{\text{d}t}={K}_{i}\left(\frac{RM}{V}\right)+{K}_{p}\left(\frac{M{\lambda}_{0}}{V}\right)-{K}_{t}\left(\frac{{\lambda}_{0}{\lambda}_{1}}{V}\right)$ (5)

$\frac{\text{d}{\lambda}_{2}}{\text{d}t}={K}_{i}\left(\frac{RM}{V}\right)+{K}_{p}M\left(\frac{{\lambda}_{0}+2{\lambda}_{1}}{V}\right)-{K}_{t}\left(\frac{{\lambda}_{0}{\lambda}_{2}}{V}\right)$ (6)

$\frac{\text{d}{\mu}_{0}}{\text{d}t}={K}_{t}\left(\frac{{\lambda}_{0}^{2}}{V}\right)$ (7)

$\frac{\text{d}{\mu}_{1}}{\text{d}t}={K}_{t}\left(\frac{{\lambda}_{0}{\lambda}_{1}}{V}\right)$ (8)

$\frac{\text{d}{\mu}_{2}}{\text{d}t}={K}_{t}\left(\frac{{\lambda}_{0}{\lambda}_{2}}{V}\right)$ (9)

where ${\lambda}_{i}$ is the momentum i of the radical species and ${\mu}_{i}$ is the momentum i of the dead polymeric species summing to;

${\lambda}_{i}=\underset{n=1}{\overset{\infty}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{n}^{i}{R}_{n}\left(mol\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mu}_{i}=\underset{n=1}{\overset{\infty}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{n}^{i}{P}_{n}\left(mol\right)$ (10)

A detailed procedure of deriving the moment equations can be found elsewhere in the works of [20] [21] .

The initial conditions used in this work: $I\left(0\right)={I}_{0}$ , $M\left(0\right)={M}_{0}$ , ${\lambda}_{i}\left(0\right)={\mu}_{i}\left(0\right)=R=0$ .

At any time and point within the system, predictions for the monomer conversion (X_{m}), number and weight average molecular weights (M_{n} and M_{w}) can be given as;

${X}_{m}=\frac{{M}_{0}-M}{M}$ (11)

${M}_{n}=M{W}_{m}\left(\frac{{\lambda}_{1}+{\mu}_{1}}{{\lambda}_{0}+{\mu}_{0}}\right)$ (12)

${M}_{w}=M{W}_{m}\left(\frac{{\lambda}_{2}+{\mu}_{2}}{{\lambda}_{1}+{\mu}_{1}}\right)$ (13)

where MW_{m} is the molecular weight of the MMA monomer. Then the polydispersity index can be written

$PDI=\frac{{M}_{w}}{{M}_{n}}$ (14)

The enthalpy change Q during the reaction progress can be defined as

$\frac{\text{d}Q}{\text{d}t}=-57700\left[-{K}_{p}\left(\frac{M{\lambda}_{0}}{V}\right)-{K}_{i}\left(\frac{RM}{V}\right)\right]$ (15)

In order to account for the cage, gel and glass effects, [12] reformulated the model of [10] basing on free volume theory and the equations below were used to define the variation of the initiator efficiency f, the termination and propagation rate constants with temperature

$\frac{1}{f}=\frac{1}{{f}_{0}}\left[1+{\theta}_{f}\left(T\right)\frac{M}{V}\frac{1}{\mathrm{exp}\left[{\epsilon}_{13}\left\{-\psi +{\psi}_{ref}\right\}\right]}\right]$ (16)

$\frac{1}{{K}_{p}}=\frac{1}{{K}_{p,0}}+{\theta}_{p}\left(T\right)\left(\frac{{\lambda}_{0}}{V}\right)\left[\frac{1}{\mathrm{exp}\left[{\epsilon}_{i3}\left\{-\psi +{\psi}_{ref}\right\}\right]}\right]$ (17)

$\frac{1}{{K}_{t}}=\frac{1}{{K}_{t,0}}+{\theta}_{t}\left(T\right){\mu}_{n}^{2}\left(\frac{{\lambda}_{0}}{V}\right)\left[\frac{1}{\mathrm{exp}\left(-\psi +{\psi}_{ref}\right)}\right]$ (18)

$\psi =\frac{\gamma \left\{\frac{{\rho}_{m}{\varphi}_{m}{\vee}_{m}^{*}}{{\epsilon}_{i3}}+{\rho}_{p}{\varphi}_{p}{\vee}_{p}^{*}\right\}}{{\rho}_{m}{\varphi}_{m}{\vee}_{m}^{*}{\vee}_{fm}+{\rho}_{p}{\varphi}_{p}{\vee}_{p}^{*}{\vee}_{fp}}$ (19)

${\psi}_{ref}=\frac{\gamma}{{\vee}_{fp}}$ (20)

In the above equations, the volume of the mixture V is determined from;

$V=M{W}_{m}\left[\frac{M}{{\rho}_{m}}+\frac{{M}_{0}-M}{{\rho}_{p}}\right]$ (21)

With
${\rho}_{m}$ and
${\rho}_{p}$ as the densities of the monomer and polymer whereas the volume of the cage, gel and glass effects of the monomer V_{fm} and the corresponding volume of the polymer V_{fp} are calculated from;

${V}_{fm}=0.149+2.9\times {10}^{-T}\left(T-273.15\right)$ (22)

${V}_{fp}=0.0194+1.3\times {10}^{-4}\left(T-273.15-{10}^{5}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}T<\left({10}^{5}+273.15\right)$ (23)

The volumetric fraction of the monomer with respect to V, denoted as ${\varphi}_{m}$ and the volumetric fraction of the polymer in respect to the same, ${\varphi}_{p}$ are obtained from

${\varphi}_{m}=\frac{\frac{M\left(M{W}_{m}\right)}{{\rho}_{m}}}{M\left(\frac{M{W}_{m}}{{\rho}_{m}}\right)+\frac{\left({M}_{0}-M\right)\left(M{W}_{m}\right)}{{\rho}_{p}}}$ (24)

${\varphi}_{p}=1-{\varphi}_{m}$ (25)

From Equations (26) to (28); ${\theta}_{f},{\theta}_{p},{\theta}_{t}$ are best fit correlation parameters that are applied as adjustable parameters for the initiator efficiency, propagation constant rate and disproportion termination constant rate respectively and may be calculated from

${\mathrm{log}}_{10}\left[{\theta}_{t}\left(T\right)\text{s}\right]=1.241\times {10}^{2}-1.0314\times {10}^{5}\left(\frac{1}{T}\right)+2.2735\times {10}^{7}\left(\frac{1}{{T}^{2}}\right)$ (26)

${\mathrm{log}}_{10}\left[{\theta}_{p}\left(T\right)\text{s}\right]=8.03\times {10}^{1}-7.50\times {10}^{4}\left(\frac{1}{T}\right)+1.765\times {10}^{7}\left(\frac{1}{{T}^{2}}\right)$ (27)

${\mathrm{log}}_{10}\left[{10}^{3}\times {\theta}_{f}\left(T\right){\text{m}}^{\text{3}}\cdot {\text{mol}}^{-1}\right]-40.8695+1.7179\times {10}^{4}\left(\frac{1}{T}\right)$ (29)

The ratios of the critical volumes of monomer i = 1 and initiator i = I over the polymer are determined as

${\epsilon}_{13}=\frac{{\vee}_{m}^{*}M{W}_{m}}{{\vee}_{p}^{*}{M}_{jp}}$ (29)

${\epsilon}_{I3}=\frac{{\vee}_{I}^{*}M{W}_{I}}{{\vee}_{p}^{*}{M}_{jp}}$ (30)

where ${\vee}_{m}^{*}$ , ${\vee}_{I}^{*}$ , ${\vee}_{p}^{*}$ are the critical volumes for the monomer, initiator and polymer respectively. The molecular weights for the monomer, polymer and initiator are denoted as $M{W}_{m}$ , ${M}_{iP}$ , $M{W}_{I}$ respectively. The parameters used in solving the above equations are given in Table 2.

Table 2. Parameters used in this study [12] [13] .

3. Results and Discussion

A program to solve the above equations was developed using the object-oriented language Python 3.5 installed via the Anaconda distro. The NumPy module [22] and a customized module were used to solve for the stiff ODEs and results were ported to matplotlib whose output has a close resemblance to the well-known MATLAB format. In order to simulate the performance of the model, it was assumed that 2 tonnes of monomer were fed to the reactor and the initiator (1% wt. of monomer) was added in 1 minute at constant rate after which the addition rate was stopped. The temperature was fixed at 70˚C. The reaction progressed for the 3 hours and the reactor profiles in Figures 1-8 were obtained.

Figure 1. Monomer conversion versus time.

Figure 2. Rate of change of initiator as the reaction progresses.

Figure 3. Monomer consumption with time.

Figure 4. Enthalpy change with time.

There is a fairly linear increase in the monomer conversion with time until the start of the gel effect as observed in Figure 1. The reaction rate is then auto-accelerated thus shifting the reaction to the limiting conversion. It should be noted that maximum monomer conversion and polymer production are limited by short reaction batch times, the existence of diffusional effects in the reaction medium and by the fast decay of the initiator as observed in Figure 2.

Figure 5. Variation of Mn and Mw with conversion.

Figure 6. Weight average molecular weight versus time.

As the initiator is added during the first minute at a constant rate, there’s an initial increase at loading time but then the initiator and monomer amounts (see Figure 3) decrease with time when the initiator radicals combine with the monomer particle leading to chain propagation which finally results into formation of the polymer.

The reaction between the monomer and initiator is highly exothermic leading

Figure 7. Number average molecular weight versus time.

Figure 8. Polydispersity index versus time.

to large amounts of heat being released. Figure 4 shows the enthalpy change as the reaction progresses. It is observed that there is an increase in the enthalpy change with time a scenario which may be attributed to the changing reaction rate constants that vary with temperature at a given time. However a maximum enthalpy change is attained prior to the depletion of the monomer inside the reactor at conversions of about 92% as shown in Figures 1-3.

For the polymer being produced, the molecular weight predictions of the simulated model (Mn, Mw and PDI) are given in Figures 5-8. It can be seen that the there is a slight drop in the polymer molecular weight due to a contraction in volume at low conversions.

However on the onset of the gel effect, larger polymer chains are produced which leads to a rise in the predicted molecular weights. As monomer conversion increases, the viscosity of reacting mixture increases. This causes a cage effect till the polymer chain stops growing. The polydispersity index predicted in Figure 8 is 27, a value greater than one indicating that the produced polymer is polydisperse.

4. Conclusion

The simulation of suspension polymerization model has been performed in Python 3.5. The prediction results were in good agreement as those reported in literature. A maximum monomer conversion of about 92.8% at a minimum batch time of about 2.2 hours was achieved at the specified conditions. The results of this work have demonstrated that Python can be employed to perform modeling and simulation studies of polymerization based processes effectively with equal success as any other programming language.

References

[1] Arora, P., Jain, R., Mathur, K., Sharma, A. and Gupta, A. (2010) Synthesis of Polymethyl Methacrylate (PMMA) by Batch Emulsion Polymerization. African Journal of Pure and Applied Chemistry, 4, 152-157.

[2] Billmeyer, F. W. (1971) Textbook of Polymer Science. John Wiley & Sons, New York.

[3] Odian, G. (1970) Principles of Polymerization. McGraw-Hill, New York.

[4] Allen, G. and Bevington, J.C. (1989) Comprehensive Polymer Science. Pergamon Press, New York.

[5] Caytano, T.R.A. and Machado, R.A.F. (2011) Soft-Sensor Based on Artificial Neuronal Network for the Prediction of Physico-chemical Variables in Suspension Polymerization Reactions. Chemical Engineering Transactions, 529-534.

[6] Verros, G.D., Achilias, D.S. and Giannoukos, G.I. (2011) Development of a Comprehensive Mathematical Model for Free Radical Suspension Polymerization of Methyl Methacrylate. Polymer Engineering and Science, 51, 670-678.
https://doi.org/10.1002/pen.21865

[7] Klaus, A., Manfred, S. and Thomas, R. (2016) Polymethacrylates. in Ullmann’s Polymers and Plastics: Products and Processes, Wiley-VCH, 885-897.

[8] Silvia, C. and Victor, B. (1999) Free Radical Polymerization of Methyl Methacrylate Modeling and Simulation under Semibatch and Nonisothermal Reactor Conditions. Journal of Applied Polymer Science, 74, 2561-2570.

[9] Chiu, W.Y., Carratt, G.M. and Soong, S.D. (1983) A Computer Model for the Gel Effect in Free-Radical Polymerization. Macromolecules, 16, 348-357.

[10] Achilias, D.S. and Kiparissides, C. (1992) Development of a General Mathematical Framework for Modelling Diffusion-Controlled Free-Radical Polymerization Reaction. Macromolecules, 25, 3739-3750. https://doi.org/10.1021/ma00040a021

[11] Garg, D.K., Serra, C.A., Hoarau, Y., Parida, D., Bouquey, M. and Muller, R. (2014) Analytical Solution of Free Radical Polymerization: Applications Implementing Gel Effect Using AK Model. Macromolecules, 47, 7370-7377.
https://doi.org/10.1021/ma501413m

[12] Seth, V. and Gupta, S.K. (1995) Free Radical Polymerizations Associated with the Trommsdorff Effect under Semibatch Reactor Conditions: An Improved Model. Journal of Polymer Engineering, 15, 283-323.
https://doi.org/10.1515/POLYENG.1995.15.3-4.283

[13] Ghosh, P., Gupta, S.K. and Saraf, D.N. (1998) An Experimental Study on Bulk and Solution Polymerization of Methyl Methacrylate with Responses to Step Changes in Temperature. Chemical Engineering Journal, 70, 25-35.
https://doi.org/10.1016/S1385-8947(98)00064-3

[14] Rosa, J.N.D. (2011) Modeling and Simulation of a Lab-Scale Polymerization Reactor. Disserta??o de mestrado, ROSA, Jaquelino Nascimento da.

[15] Ekpo, E.E. (2006) Dynamic Optimisation and Control of Batch Polymerization Process. Ph.D. Thesis, University of Bradford, Bradford.

[16] Fangbin, Z., Gupta, S.K. and Ray, A.K. (2001) Modeling of the Sheet-Molding Process for Poly (Methyl Methacrylate). Journal of Applied Polymer Science, 81, 1951-1971.

[17] Moran, S. (2015) An Applied Guide to Process Plant Design. Butterworth-Heinemann, Oxford.

[18] Achilias, D.S. and Kiparissides, C. (1994) On the Validity of the Steady State Approximations in High Conversion Diffusion Controlled Free-Radical Copolymerization Reactions. Polymer, 35, 1714-1721.
https://doi.org/10.1016/0032-3861(94)90846-X

[19] Penlidis, A., Ponnuswamy, S.R., Kiparissides, C. and O'Driscoll, K.F. (1992) Polymer Reaction Engineering: Modeling Consideration for Control Studies. The Chemical Engineering Journal, 50, 95-107.
https://doi.org/10.1016/0300-9467(92)80013-Z

[20] Dotson, N.A., Galvan, R., Laurence, R.L. and Tirrell, M. (1996) Polymerization Process Modeling. VCH Publishers, New York.

[21] Ray, W.H. (1972) On the Mathematical Modeling of Polymerization Reactors. Journal of Macromolecular Science, Part C, 8, 1-56.

[22] Kiusalaas, J. (2010) Numerical Methods in Engineering Using Python. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511812224