Stochastic Simulation of Droplet Interactions in Suspension Polymerization of Vinyl Chloride

In this paper a population balance based mathematical model is presented for describing suspension polymerization of vinyl chloride. The properties of the polymer product and the behaviour of the stirred batch polymerization reactor are investigated by simulation. Two-phase kinetics model of free radical polymerization is used, and heat balance is also included into the model. Beside the coalescence and breakage phenomena, are taken interchanges of species and heat between the droplets induced by collisions into account forming a complex threescale system. The motion of droplets in the physical space of the polymerization reactor, the breakage, coalescence and coalescence/redispersion processes are simulated by using a coupled continuous time – Monte Carlo method.


Introduction
Suspension polymerization is widely used for commercial manufacture of many important polymers including poly(vinyl chloride) (PVC), poly(methyl methacrylate), expandable polystyrene, styrene-acrylonitrile copolymers and a variety of ion exchange resins.
Suspension polymerization is usually carried out in batch tank reactors.In suspension polymerization, one or more water-insoluble monomers containing oil-soluble initiator(s) are dispersed in a continuous liquid phase by a combination of strong stirring and the use of small amounts of suspending agents (stabilizers).The stabilizers hinder first the coalescence of the monomer droplets, and later stabilize the polymer beads whose tendency to agglomerate may become critical when the polymerization has advanced to the point where the polymer beads become sticky [1].PVC is produced by powder polymerization since PVC is insoluble in the vinyl chloride monomer (VCM) hence it immediately precipitates out formatting a separate phase.In suspension polymerization of VCM the free radical polymerization reactions take place in the monomer droplets because the initiator is dissolved in those.
The properties of PVC are influenced by the polymerization conditions: polymerization temperature, stirring conditions, reactor size, concentration of initiator, adding type and number of initiators.The temperature of polymerization is one of the most important parameters.Only 1 or 2 K differences in temperature can induce significant changes in the properties of polymer products.
Two cases can be distinguished as regards the addition of initiator.In the first case, the monomer and initiator are homogenized first, the mixture is filled into the reactor, and subsequently the reactor is heated to the reaction temperature.The advantage of this process is that the distribution of initiator in the monomer is perfect while its disadvantage is that the polymerization reactions occur under non-isotherm conditions during the heating stage.Therefore, the properties of the end-use polymer might be not uniform.In the other case, first the monomer and water are filled into the reactor, the suspension is heated to the polymerization temperature and thereafter the initiator is added to the suspension.The advantage of this method is that the polymerization reactions occur under isothermal conditions and the properties of the end-use polymer usually are uniform.The disadvantage of this process is that the initiator distribution in the monomer droplets may be not perfect.In this work we analyse the second case.
Two processes occur parallel with each other inside the suspension polymerization reactor.On the one hand the highly exothermic polymerization reactions inside the monomer droplets, having rates depending on the actual reactants concentrations and droplets temperature, form a continuous in time deterministic process.Simultaneously with that the meso-scale interactions take place in the reactor form discrete event processes: random droplets coalescence, binary breakage of droplets and binary collisions (coalescence/redispersion) may occur between the droplets moving in the reactor space the frequencies of which depends on the size and number of droplets in a unit volume of the reactor.
The physical properties in a suspension polymerization reactor significantly affect the droplet size distribution.Droplet size depends on the physical properties of the two phases, the phase ratio, the nature of the suspension flow, and the condition of the phase interface.Droplet breakage in agitated suspensions can be caused either by frictional forces (through viscous shear) or by inertial forces (through turbulence) [2].In industrial suspension polymerization reactors, the volume fraction of the dispersed phase is usually high and droplet breakup is accompanied by droplet coalescence.Thus, the average droplet size and the droplet size distribution are both influenced by droplet breakage and droplet coalescence.Thus, the framework of population balances is ideally suited to the description of the complex dynamics of particulate polymerization processes [3].
Population balance models for modelling the suspension polymerization were applied by Alexopoulos and Kiparissides [4], Kiparissides [5] and Alvarez et al. [6].All these models have taken only the droplet size distributions and its variation into account.However, in suspension polymerization, because of changes of droplet volume, concentrations of species and even of the temperature of droplets are important thus the population balance equation becomes, in principle, multi-variable.
Solution of multi-variable population balance equations is a crucial problem.As analytical solutions of PBEs are available in very few cases, numerical techniques are essential in most practical applications.There are several numerical methods available that meet the accuracy requirements [7] but, taking the random nature of breakage, coalescence and micromixing processes into consideration Monte Carlo (MC) simulation also seems to be a suitable method for solving this numerical problem.
However, in the case of suspension polymerization besides the multi-variable problem there exists another problem.Namely, here time continuous processes and discrete event processes occur in parallel thus solution of the system of equations requires a special technique -combination of a continuous time and discrete event treatment.This phenomenon was studied by Salikas et al. [8] but they didn't describe the problem explicitly.
Irizzary [9,10] also presented a case in which there were chemical reactions in the dispersed phase.However, Irizarry assumed instantaneous chemical reaction in the droplets while in polymerization reactions take hours but are highly exothermic thus in droplets should be taken the heat effects into consideration.
In a model of suspension polymerization reactor should be taken both processes, i.e. the deterministic polymerization reactions and the discrete event processes of meso-scale interactions into consideration and this can be achieved by using the population balance approach.
During the course of polymerization the volume of droplets is also changed due to the density difference of monomer and polymer and due to the coalescence and breakage of droplets that should be taken into account as well.In suspension polymerization of vinyl chloride, however, the volume fraction of the dispersed phase is about 0.3 therefore the total volume of mixture in suspension polymerization can be considered constant.
The aim of the present paper is to model and analyse a suspension polymerization process of vinyl chloride in which beside the polymerization reactions inside the droplets collision induced heat and mass exchange, coalescence of droplets and breakage of droplets occurs in the reactor.The kinetic data of vinyl chloride polymerization are taken from the literature [11].The effect of initial droplet size distribution, the effect of initial initiator distribution and the effect of droplet interactions are analysed under isothermal and non-isothermal conditions on monomer conversion and on some average polymer properties.

Multi-dimensional population balance model
Consider a large population of interacting monomer droplets moving stochastically in the continuous carrier phase of an intensively stirred reactor.The initiator, soluble in the monomer is distributed in the droplets also by stirring.As the reactor is heated to the polymerization temperature the polymerization process is started under isothermal conditions and subsequently two processes occur parallel in the reactor.Polymerization reactions [12], having rates depending on the actual states of phases and concentrations inside the droplets form a continuous in time deterministic process, while coalescence and breakage of droplets with micromixing in the coalescence state form a stochastic discrete event process.
Let υ denote the volume coordinate of droplets, c stand for the vector of concentrations of K ≥ 0 relevant chemical species inside the droplets, and T denote the temperature of droplets.When analysing the properties of the polymer product beside tracking the changes of concentrations of the initiator and monomer it is reasonable to compute also the first three leading moments of the live and dead polymer chains.This formulation requires 8 variables and 8 differential equations and provides a sufficiently detailed description of the polymerization reactions [12,13,14].Therefore in suspension polymerization of vinyl chloride K=8 and the vector of concentrations c = (c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , c 7 , c 8 ) denotes, in turn, the concentrations of the initiator and monomer, and the three leading moments of the live and dead polymer chains in a droplet.
Then, assuming that the reactor is perfectly mixed at macroscale and the motion of droplets in the continuous carrier phase is fully stochastic without orientation the micro-scale state of a droplet, in general case, is given by the vector (υ,cT)  R K+2 without indicating its position and velocity in the space of the reactor [15].Defining the population density function given as a mapping (υ, c, T, t) → n(υ, c, T, t) by means of which n(υ, c, T, t)dυdcdT provides the number of droplets from the volume (υ, υ+dυ) and temperature (T, T+dT) intervals and concentration region (c, c+dc) at time t in a unit volume of the suspension.
In suspension polymerization the concentrations of spesies in the droplets change continuously in time because of the polymerization and may also change jump-like in both pure mass exchange and aggregation events are collision induced.The concentration of droplets in a break up event is assumed to be unchanged.The volume of droplets is also changed continuously in time because of the polymerization reactions and the significant difference between the densities of the monomer and the polymer but their volumes may also change jump-like because of collision-induced aggregation and breakage.The temperature of droplets is also changed continuously in time because of the polymerization reactions but their temperatures may also change jump-like because of collision-induced events, the pure heat exchange and aggregation events.The temperature of droplets in a break up event is assumed to be unchanged.
The meso-scale model of the droplet population in the polymerization reactor takes the form where the rates of change of the population density function because of the deterministic continuous processes are on the left hand side while the terms on the right hand side represent the random collision-induced processes such as collision induced mass and heat exchange between the droplets, termed coalescence/redispersion, droplets coalescence and breakage.These terms were detailed and analysed by Lakatos [15], we present the solution of the population balance equation (1) here using a Monte Carlo method coupled with the continuous time treatment of deterministic processes.

Solution by means of a coupled continuous time -Monte Carlo method
MC methods can be divided into two classes according to the treatment of the time step.These are referred to as "timedriven" and "event-driven" MC [10].Here we used an eventdriven MC method as follows.Breakage and coalescence of droplets are affected by the droplet volumes and the local mechanical conditions in the dispersion, turbulent energy dissipation and shear forces as well as by physical properties: viscosity and density of the phases, interfacial tension, and other interfacial phenomena, such as the surface charge of droplets.The simulated droplet diameter interval is sectioned to equal parts termed size classes.All droplets are assigned into these classes based on their actual diameters.The volumes of droplets are changed continuously in time because of the polymerization reactions and the significant difference between the densities of monomer and polymer and their volumes can be changed jump-like because of collision-induced aggregation, coalescence/redispersion and breakage.Because the change of volume of droplets the number of droplets in classes can also vary.The collision and breakage processes are modelled as inhomogeneous Poisson processes independent from each other the intensities of which are computed individually for all classes.Naturally, in the time intervals between subsequent events computations are carried out using homogeneous approximations.The droplet or droplets which take a share in events are selected randomly using random numbers generated from the uniform probability distribution in (0, 1).The steps of the solution method are as follows: Initialization: Initial droplet size distribution is generated using the beta distribution and all state variables are given initial value.The number of droplets is N. Set the time equal to zero.The investigated droplet diameter interval is 1 µm to 250 µm and this interval is sectioned to 25 equal parts.All state variables are given initial values (initiator concentration, temperature of polymerization, etc.).
Step 1: Selection of the next event and next event time from all possible events using equations and where λ k denote the intensities (mean frequencies) of the collision and break up events computed individually for all size classes.Set the simulation time to: (1) Step 2: The polymerization reactions take place in the monomer droplets, so, for all droplets integrate the set of intraparticle reactions from t i to t i+1 .The concentration, temperature and volume of droplets are changed continuously in time because of the reactions.
Step 3: The formerly selected event occurs.The event is collision of two droplets.aa) If this event is coalescence/redispersion then two droplets are selected randomly from the diameter classes d and d'.A random number ω k  [0,1] is generated to calculate the rate of species mass and heat exchange between the colliding droplets.The total number of droplets does not change and remains N. ab) If this event is coalescence then two droplets are selected randomly from the diameter classes d and d'.These droplets are eliminated and a new droplet is formed from these droplets with size υ new = υ i + υ j .The properties i.e. concentrations of initiator and monomer, conversion, moments and temperature of the new droplets are calculated from the properties of coalesced droplets assuming homogeneous distributions of intensives.Then set N:=N-1.The event is break up of a droplet.A single droplet is selected randomly from the diameter class d.This droplet is eliminated and two new twin droplets are formed from that with υ new = υ i / 2 .All extensive quantities are transformed in correspondence with the volume while the intensives become homogeneous of the same values as that of the mother one.Then set N: = N + 1.
Step 4: If t i+1 ≥ t final then stop, i.e. end the simulation, otherwise go to Step 1.

Results and discussion
The corresponding computer program and all simulation runs were written and carried out in MATLAB environment.
The most commonly used model developed for modelling the suspension polymerization of VC is the two-phase model [11,16].
The key feature in all these models is that PVC is practically insoluble in its monomer, and polymerization proceeds simultaneously in the two phases almost from the start of the reaction.These models assume from kinetic point of view the polymerization of VC is considered to take place in three stages [17].
Stage 1.During the first stage the droplets contain as good as completely pure monomer, although PVC is insoluble in its monomer, but under 0.1% monomer conversion the polymer concentration is below the solubility limit.In this stage we calculate the mass balance equations only for the monomer rich phase.
Stage 2. This stage begins from the appearance of the separate polymer phase to the so called critical conversion, X c , at which the separate monomer phase disappears (0.1% < X < X c ).The free radical polymerization reactions take place in both the monomer and polymer rich phases at different rates and is accompanied by transfer of monomer from the monomer phase to the polymer phase so that the latter is kept saturated with monomer.The disappearance of the monomer phase is associated with a pressure drop in the reactor.During this stage we calculate the equations for both of phases.It means that the number of equations is doubled.The reaction rate coefficients are different in the two phases.
Stage 3. Above the critical conversion the polymerization reactions take place only in the polymer rich phase.The monomer mass fraction in the polymer phase decreases as the total

Yes
No monomer conversion approaches a final limiting value.In this stage we calculate the balance equations only for the polymer rich phase.Since polymerization is exothermic, therefore the heat balance in droplets was calculated, too.The energy balance for a droplet is given with the following equations: where and In Eqs (4-6), ρ, c p , ΔH r , m 0 , M w , X, h, A, T d and T c denote, in turn, the density of droplets, their heat capacity, total enthalpy change of reactions, initial mass of monomer, molecular weight of monomer, monomer conversion, the heat transfer coefficient and heat transfer surface area of droplets emerged in the continuous phase, the temperature of droplets and temperature of the continuous phase.
In all simulation runs, it was assumed that the temperature of the continuous phase, because of a well-controlled cooling process was constant while the heat transport between the droplets and the continuous phase was computed using the correlation by Ranz and Marshall [18].Some preliminary simulations were carried out in which the effect of the total number of elements was analysed.These simulation runs revealed that a 1000-5000 elements droplet population proved to be sufficient to simulate the process with reliable approximation [14].
In the first step the simulation program was identified.The kinetic data and model equations of polymerization kinetic were taken from literature [11].In this case, the droplets were mono-dispersed as regards their volume while the initial distribution of the equal amounts of the initiator in the droplets was uniform.As a consequence, this case provided an ideal reactor mixed perfectly both at macro and micro levels, and all droplets behaved in the system as identical perfectly mixed micro-reactors.Figure 2 shows the temporal evolution of monomer conversion.It is seen that the time profile obtained by simulation fit very well to the experimental.So this computer program proved to be suitable for simulating the suspension polymerization of VC.
The effect of initial initiator distribution under isothermal conditions was analysed in our previous works [12,13,14].In these papers the initial droplet size distribution was uniform, and we assumed that the temperature of the monomer droplets and the continuous phase were constant which means the cooling of the reactor was perfect.In our other work we analysed the effect of the temperature rise of droplets [19].In this study the initial droplet size distribution was uniform, the computations was carried out under conditions that all droplets contain initiator but not the same amount.Assuming perfect stabilization, coalescence and breakage of droplets becomes negligible.Besides, taking into account that the rate of diffusional mass transfer is smaller with some orders of magnitude than that of heat it is justified to assume that collisions of droplets, termed collision/redispersion events hereafter, induce only heat exchange between the colliding droplets, caused by their possible temperature differences, with negligible species mass interactions.Afterwards we analysed the effect of initial initiator distribution generating initial droplet size distribution at the beginning of simulation [20] and along with the calculation of temperature rise of droplets [21].In these simulations were taken into consideration all of meso-scale interactions, such as coalescence, breakage and coalescence/redispersion of droplets.
In our present work we analysed the effect of initial droplet size distribution and the effect of meso-scale interactions during polymerization.In the first case the initial droplet size distribution was uniform which means all droplets had the same size and volume.In the other case the monomer droplets had initial droplet size distribution.The initial droplet size distribution is affected by the collision, coalescence and breakage frequencies.These frequencies were calculated using the expressions and parameters by Alopaeus et al. [22] while the energy dissipation rate, influenced by some properties of the impeller such as impeller speed, impeller diameter, power number of impeller, and naturally the properties of mixed phase, the volume and density of the mixture was computed as discussed by Nere et al. [23].
Figure 3 represents the monomer droplet size distribution computed for impeller speed is 350 rpm and dispersed phase volume fraction is 0.3.Comparing that with the monomer droplet size distribution measured by Zerfa and Brooks [24] for the same parameter values the simulation data fit in really well with the measured data as it is seen in Fig. 3. Experimental data (Sidiropoulou and Kiparissides,1990) Calculated data Fig. 2 The time profiles of experimental data from literature and simulation data.The droplet size distribution was uniform, all droplets contained the same amount of initiator, the amount of initiator was 0.0029 mole % based on monomer and the temperature was 323.
First we analysed the effect of initial droplet size distribution at different initial initiator distribution.The temperature of the continuous phase was constant, 323 K.We compared the cases when initially all droplets have the same size to the case when the droplets have initial monomer droplet size distribution (Fig. 3b).In the case when all droplets have the same size the droplet diameter was 50 µm.During these investigations were taken all meso-scale interactions, the coalescence, the binary breakage and the coalescence/redispersion of droplets into consideration.Figure 4 shows the effect of initial monomer droplet size distribution at different initial initiator distributions.Namely, initial the random 25 % (Fig. 4a), 75 % (Fig. 4b) or 100 % (Fig. 4c) of droplets contain the same amount of initiator.
From Fig. 4 we can see that the effect of initial monomer droplet size distribution is significant in the case, when initially only less than 75% of droplets contain initiator.In Fig. 4a we can see that the difference between the cases is significant.Figure 5 represents the number average molecular weights in case when initially 25% of droplets contain initiator.We can see that the time profiles of the two cases are really similar, but in the middle of the process there is a small difference in the profiles.
The simulation results (Fig. 4 and 5) demonstrated that the ideal case is when initially the monomer droplet size distribution is uniform.
We analysed the effect of the droplet interactions, the temperature of the continuous phase was 323 K.During the simulations the monomer droplets had initial droplet size distribution (Fig. 3b).We compared two cases, in the first one were taken all meso-scale interactions into consideration, while in the second case we assumed that the stabilization of droplets was perfect, therefore the coalescence and binary breakage of droplets are negligible and only coalescence/redispersion of droplets takes place in the reactor.
From Fig. 6 we can see that the effect of meso-scale interactions is significant in case, when initially only less than 75 % of droplets contain initiator.In Fig. 6a and Fig. 6b we can see that differences between the cases are significant.Figure 7 represents the number average molecular weights in cases when initially 25 % (Fig. 7a) and 75 % (Fig. 7b) of droplets contain initiator.
We can see that the time profiles in the case when initially 75% of droplets contain initiator are really similar, but at the middle of the process there is a small difference in the profiles.In the case when initially only 25% of droplets contain initiator the difference between the two cases is remarkable.The end use polymer properties are different in this case.
The simulation results (Fig. 6 and 7) demonstrated if the droplets initially have droplet size distribution the ideal case is when all events (coalescence, coalescence/redispersion, breakage) take place in the reactor.It is because all events can help to distribute the initiator between droplets.

Conclusions
A multi-variable population balance model was developed for suspension polymerization of vinyl chloride in a batch polymerization reactor.The model contain two parallel processes, the polymerization reactions inside the monomer droplets form a continuous in time deterministic process and the meso-scale interactions of droplets form discrete event processes.The heat balance of monomer droplets was taken into the model.A coupled continuous time -Monte Carlo method was developed for the solution of a multi variable population balance equation.A "time-driven" MC method was coupled with the deterministic continuous time computations of polymerization reactions.The breakage and collision frequencies were calculated according to the diameter classes of monomer droplets and the number of droplets inside the classes.The next event and next event time was calculated from these frequencies and the polymerization reactions were calculated during this time as deterministic continuous time processes.
The effect of initial droplet size distribution and the effect of meso-scale interactions were analysed on monomer conversion and on some average polymer properties.The simulation results demonstrated if the initial droplet size distribution is uniform then the monomer conversion is higher at same initial initiator distribution.If the initial droplet size distribution is not uniform the ideal case is when all meso-scale interactions take place in the reactor, because all these events help to distribute the initiator between the droplets."National Excellence Program -Elaborating and operating an inland student and researcher personal support system".This project was subsidized by the European Union and co-financed by the European Social Fund.The financially supports are gratefully acknowledged.initial droplet size distribution, all events initial droplet size distribution, c/r a., Initially 25% of droplets contain the all amount of initiator b., Initially 75% of droplets contain the all amount of initiator Fig. 7 The effect of meso-scale interactions on number average molecular weight, initially 25% of droplets contain initiator.The amount of initiator was 0.0029 mole % based on monomer and the continuous phase temperature was 323 K

Fig. 1
Fig. 1 Algorithm of coupled continuous time -Monte Carlo method

Fig. 3 Fig. 4
Fig.3 Comparison of experimental[21] and calculated data of the initial droplet size distributions.The impeller speed was 350 rpm, the dispersed phase volume fraction was 0.3.

Fig. 6 Fig. 5
Fig.6 The effect of meso-scale interactions on monomer conversion at different initial initiator distribution.The amount of initiator was 0.0029 mole % based on monomer and the continuous phase temperature was 323 K