Temperature profiles in batch methyl methacrylate polymerization in gelled suspension.
Suspension polymerization is a well-known process for the production of small polymeric particles. It offers two main advantages: it avoids thermal problems by providing a suspending phase, which maintains a low viscosity also at high conversion, permits complete agitation during the whole process, and absorbs part of the heat produced. Moreover, it produces particles of homogenous size and shape, which are ready to use in different applications; these advantages are only attained if agitation is maintained during polymerization and a suspending agent is added, which guarantees the absence of coalescence when reaction proceeds. Therefore, if the polymer is to be used for biomedical purposes, expensive purification from the suspending agent is required.
This is why a process has recently been developed for methyl methacrylate (MMA) polymerization (1, 2) to produce polymeric particles, e.g., to be used as bone cement. In this process the particles are of the same size and shape as in the usual suspension process, but no suspending agent is required. This is attained by a polymerization into small particles embedded in an unstirred, gelled, suspending phase. The monomeric particles are formed by stirring a mixture of water, agarose, monomer, and initiator at a temperature at which polymerization is prevented by the slow kinetics and gelification is prevented by the vigorous stirring. Once the particles have been formed, the operator stops the mixer, and the suspending phase, composed of water and agarose, gels. The jacket temperature at this point is raised to a level at which polymerization can start in the imprisoned droplets, while no coalescence can occur in the gelled phase. In this way only one of the main advantages of the suspension process (the one related to the size and shape of the particles) is maintained while the heat dissipation is partially prevented by the lack of agitation. The process feasibility is therefore related to the thermal behavior of the reactor and in particular to the maximum temperature peak.
A feasibility study by fundamental modeling has been recently developed (3); however, this work, mainly interested in validating the kinetic and physical model, requires support by a sensitivity analysis. This is meant to test if the reactor operates at or near conditions for which the risk of runaway is not present and if small errors in the identified physical and kinetic parameters are not largely amplified in the system state predictions.
Parametric sensitivity has inspired extensive research and review (4-6). It is likely that the first work in this field was by Semenov (7), regarding thermal explosion theory. In the following years, different techniques were developed, based mainly on intuitive criteria. A practical approach [e.g., Rawlings and Ray (8)] consists in simply changing the parameters one at a time by a reasonable amount and examining the system behavior in the modified conditions. In this way, variations can be analyzed and the designer can decide which operating variables can be changed to obtain the desired product, which kinetic and physical parameters can be chosen to improve the accuracy of the system description, and in which conditions the system should not operate for high sensitivity problems.
More rigorous studies (9-12) eliminate the arbitrary choice of how much the parameters should be changed to study the modifications in the system states.
All these methods are related in different ways to the partial derivatives of the states of interest with respect to the system parameters.
Morbidelli and Varma (9) proposed a technique that they applied to thermal explosion theory for the computation of the sensitivity of the value of conversion that corresponds to the maximum temperature. Their method does not add much complexity to the mathematical description of the system; moreover, once a parameter has been taken into account, the computation of sensitivity for other parameters is straightforward. The same technique has been applied in later works to chemical processes, including polymerization with (10) and without (11) a gel effect. A different approach is the one of Caracotsios and Stewart (12), who propose a generalized technique for initial value problems described by systems of coupled differential and algebraic equations; sensitivity equations are directly solved via the local Jacobian of the state equations. Application of the technique to a stiff industrial reactor model is presented by the same authors. Their approach is followed in this work. Software available from the authors is used both to solve the coupled partial differential equations that describe the system dynamics and to compute the sensitivities with respect to the system parameters.
The main interest in the system states concerns the temperature profiles. In particular, temperature must not exceed, anywhere in the reactor, the maximum value at which stability of the gel is guaranteed ([approximately]90 [degrees] C). Moreover, temperature is a sensitive indicator of any possible runaway in the reactor. The importance of temperature profiles in the polymerization of MMA, related to autoacceleration, has been pointed out (13).
The parameters, with respect to which the temperature sensitivities are computed, are divided into three different classes: operating variables, kinetic parameters, and physical parameters. The first concern is the choice of the initial conditions and of the jacket temperature. The kinetic parameters are related to the steps of the polymerization and to the gel and glass effects that take place as conversion increases. The last ones are related to the transport phenomena, mainly for the heat transfer.
The last two classes are of interest because it is important to understand which parameters can be optimized in order to make the model follow the process dynamics more closely, and to test if small errors in parameter estimation can cause big variations in the estimated process behavior.
On the other hand, sensitivity to the operating conditions is important both for choosing these variables properly and for understanding whether any problems in the thermal behavior of the reactor can be anticipated.
A mathematical model is described for a batch reactor for gelled suspension polymerization. It is a glass hemispheric reactor with a refrigerated jacket and a mixer. Four thermocouples at different radial positions are available to follow the system thermal behavior on-line. Figure 1 shows a schematic representation of the reactor; the gelled particles are depicted in the enlargement.
Material and Thermal Balances
In the first step of the process, suspended MMA monomeric particles are trapped into a gelled aqueous phase. The suspension is obtained in a few minutes by stirring the mixture at [approximately]40 [degrees] C, at which the decomposition of the initiator (benzoyl peroxide - BPO - has been used) is almost negligible.
The polymerization starts as soon as heat is supplied to the reacting system by means of an increase in the jacket temperature (a step variation is usual). In the outer region of the reactor, the temperature is quickly influenced by the variation of the jacket temperature. On the other side, in the inner part of the reactor, the initial increase in temperature is much slower. The temperature wave moves then gradually from the outside to the inside of the reactor and makes the reaction start. Every particle therefore undergoes a different thermal history, exchanging heat with the surrounding environment. A single particle, however, is considered a closed system from the material point of view; therefore it can be modeled as a microbatch homogeneous reactor for bulk polymerization. Since no material flux between the phases takes place, it is not considered in the model, and the possibility of reaction-diffusion instabilities is eliminated.
The kinetic scheme considered is the classical one for free radical polymerization:
[Mathematical Expression Omitted]
[Mathematical Expression Omitted] Initiation
[Mathematical Expression Omitted]
[Mathematical Expression Omitted]
Termination by disproportionation
[Mathematical Expression Omitted]
Termination by combination
where I is the initiator, R[prime] is the primary radical, M is the monomer, [P[prime].sub.n] is the growing radical of length n, [P.sub.n] is the corresponding dead polymer, and [k.sub.d], [k.sub.i], [k.sub.p], [k.sub.td], [k.sub.tc] are the kinetic constants. Chain transfer to monomer is negligible for MMA, and no other chain transfer agent is present.
Material balances for the reacting species can be written following classical models. Instead of writing individual balances for the growing radicals and the dead polymers of all lengths, the moments of the molecular weight distribution are introduced and global balances are written for them. The moment of order i of the growing chains and of the dead polymers are respectively defined as follows:
[[Lambda].sub.i] = [summation of] [n.sup.i][P[prime].sub.n] where n=1 to [infinity] [[Mu].sub.i] = [summation of] [n.sup.i][P.sub.n] where n=1 to [infinity] (1)
Balances for the monomer, the initiator, and the first two moments of the growing radicals take the following form:
[Delta]I/[Delta]t = -[k.sub.d]I - [Epsilon]I/1 + [Epsilon]x [[Lambda].sub.0](1 - x)[k.sub.p] (2)
[Delta]x/[Delta]t = [k.sub.p](1 - x)[[Lambda].sub.0] (3)
[Mathematical Expression Omitted] (4)
[Delta][[Lambda].sub.1]/[Delta]t = -[Epsilon][[Lambda].sub.1][[Lambda].sub.0]/1 + [Epsilon]x (1 - x)[k.sub.p] + 2f[k.sub.d]I - [k.sub.t][[Lambda].sub.0][[Lambda].sub.1] + [k.sub.p][[Lambda].sub.0][M.sub.0] 1 - x/1 + [Epsilon]x (5)
where [k.sub.t] = [k.sub.tc] + [k.sub.td], [M.sub.0] is the initial monomer concentration, x is the monomer conversion, f is the initiator efficiency, [Epsilon] = ([[Rho].sub.m] - [[Rho].sub.p])/[[Rho].sub.p] is the volume expansion factor, and [[Rho].sub.m] and [[Rho].sub.p] are the monomer and polymer densities, respectively. Balances for the moments of the dead polymers and for the higher moments of the growing radicals can be written analogously [see, e.g., Chiu et al. (15)].
It has been shown (14) that the termination rate constant cannot be considered independent of the length of the terminating radicals for the polymerization of methyl methacrylate. The full form of the bimolecular termination rate constant [k.sub.t](i, j) cannot be found via molecular weight distribution measurements. In any case, different [k.sub.t] averages, which depend on the order of the moment, need to be included in the polymer moments balances. However, the complete dependence of the [k.sub.t] averages on temperature and conversion is not available, so they cannot be used in this study. The dependence of [k.sub.t] on conversion and temperature is the one of Chiu et al. (15) (see the following subsection for the details), which is fit for the evaluation of the zeroeth order moment but gives errors of magnitude that increase with the order of the moment (14). In any case, given that only the first two moments of the growing radical distributions are necessary to obtain the global distributions as explained below, no large errors in the distribution are expected.
Partial derivatives have been introduced in the material balances, given that all concentrations are functions not only of time but also of the radial position at which the particle to which they refer is located.
Material balances are coupled through the temperature dependence of the kinetic constants to the thermal balance. This balance is written under two main assumptions:
1) the particles are at thermal equilibrium with the surrounding gel; this quasi-steady-state assumption (QSSA) is justified by the small size of the particles and by the correspondent very high surface-to-volume ratio, which permits a rapid local heat exchange;
2) the heat flux on the top of the reactor is negligible.
The first hypothesis permits one to use a single variable T(r, t) to describe the temperature both of the reacting phase and of the gel and to describe the system thermal behavior by a global energy balance that considers the two phases as a homogeneous medium from a thermal point of view.
The second hypothesis permits one to write the energy balance in spherical coordinates, which simplifies its formulation. Such a balance therefore takes the following form:
[Rho][c.sub.p] [Delta]T/[Delta]t = 1/[r.sup.2] [Delta]/[Delta]r ([r.sup.2]K[Delta]T/[Delta]r) + [Delta][H.sub.r][[Epsilon].sub.[Sigma]][M.sub.0][k.sub.p](1 - x(r, t))[[Lambda].sub.0](r, t) (6)
where the terms on the right-hand side represent the diffusion term and the heat generation by reaction term, respectively. [Delta][H.sub.r] is the heat of reaction, [[Epsilon].sub.[Sigma]] is the initial volume fraction of the reacting phase:
[[Epsilon].sub.[Sigma]] = [([V.sub.MMA]/[V.sub.MMA] + [V.sub.[H.sub.2]O]).sub.t=0] (7)
K and [Rho][c.sub.p] are average values of the heat conductivity and of the heat capacity, respectively. The average weight heat capacity, which is considered independent on temperature, takes the form:
[Rho][c.sub.p] = [([([Rho][c.sub.p]).sub.[H.sub.2]O][W.sub.[H.sub.2]O] + [([Rho][c.sub.p]).sub.MMA] [W.sub.MMA]/[W.sub.[H.sub.2]O] + [W.sub.MMA]).sub.t=0] (8)
It is not as straightforward to give an average value for the heat conductivity. The main uncertainty is related to the physical state of the mixture; even if from a macroscopic view, no movements are detected in the gelled phase, this is not completely true on a microscopic scale, particularly when temperature increases. This has an impact on the conductivity, which increases with temperature.
At first a volume average value of conductivity is computed as:
[Mathematical Expression Omitted] (9)
where [Gamma] = [[Epsilon].sub.[Sigma]](1 + [Epsilon]x) and [K.sub.[H.sub.2]O], [K.sub.MMA], [K.sub.PMMA] are taken from the literature (16). Thereafter, a quadratic dependence has been introduced in order to make conductivity increase with temperature:
[Mathematical Expression Omitted] (10)
where [T.sub.O] is the gelation temperature and the values of [[Alpha].sub.1], [[Alpha].sub.2] that give the best fit of experimental data have been chosen. These values, for an agarose weight fraction of 0.01 in the suspending phase, are reported in Table 1 with all the other values of the kinetic and physical parameters that have been used.
[TABULAR DATA FOR TABLE 1 OMITTED]
Finally, the model consists of a system of five partial differential equations (Eqs 2-6) with the following initial and boundary conditions:
[Mathematical Expression Omitted] (11)
[Mathematical Expression Omitted] (12)
where U is the heat transfer coefficient of the reactor jacket and [T.sub.j] is the jacket temperature. The initial conditions refer to the initial state of the reactor while the boundary ones are related to the spherical symmetry of the reactor and to the heat transfer to the jacket.
Gel and Glass Effects
The gel and the glass effects are typical phenomena that appear during the course of many polymerization reactions.
The gel effect (17-19) usually appears for monomer conversion of 20% to 40% and is mainly related to the increase in viscosity of the reacting phase. The mobility of the growing chains is strongly prevented from the high viscosity so that the termination rate constant, which is diffusion controlled, drops drastically. On the other side in the first stage, the propagation rate constant does not change significantly because the small monomer molecules maintain their mobility. The increase in the ratio [k.sub.p]/[k.sub.t] causes a fast increase in the reaction rate, which makes both viscosity and the produced heat increase further. The average molecular weight, which follows the course of reaction, moves to higher values, and if the heat produced is not removed from the system, there is the risk of thermal runaway.
When high conversions are reached, the mobility of the monomer molecules decreases as well. Therefore also the propagation step (and transfer if present) becomes diffusion controlled so that the propagation rate decreases drastically. This phenomenon, known as the glass effect, slows down the polymerization rate, which tends to zero at very high conversions. The polymer produced, when the glass effect becomes significant, is characterized by a low molecular weight.
A large number of models have been proposed in the last two or so decades to describe the gel and the glass effects. The first (20-23) were mainly based on empirical relations. A deeper understanding of the physics of the phenomena, together with the availability of more powerful computing tools, has caused the development of more complex models, most of which have a semi-empirical nature (24-26). Reviews of these models are available (27-30).
Most of the developed models that describe the gel effect introduce a diffusional restriction to the termination rate constant after a characteristic break point. This leads to the division of the operating range into three different regions (before the gel effect, after the onset of the gel effect, after the onset of the glass effect), which are described by different equations. This description, which is reasonable from a mathematical point of view, is somewhat artificial from a physical point of view.
In this work, the model by Chiu et al. (15) is adopted. This model is related to a physical description of the diffusional events based on a molecular viewpoint. It has been clearly shown (31) that this model provides a good description of conversion and molecular weight both for MMA and styrene.
A basic advantage of this model is the continuity of the mathematical relationships that make it possible to use them throughout the course of reaction. This is the main reason why the model by Chiu et al. turns out to be particularly appropriate for mathematical modeling and for engineering studies such as parametric sensitivity, control system design, and parameter optimization. In recent years many authors have adopted and studied this semi-empirical model (32-35).
The constitutive equations for the propagation and termination rate constant that are used in the model are the following:
1/[k.sub.t] = 1/[k.sub.t0] + [[Theta].sub.t](T, [I.sub.0]) [[Lambda].sub.0]/exp[2.303[[Phi].sub.m]/(A(T) + B[[Phi].sub.m])] (13)
1/[k.sub.p] = 1/[k.sub.p0] + [[Theta].sub.p](T) [[Lambda].sub.0]/exp[2.303[[Phi].sub.m]/(A(T) + B[[Phi].sub.m])] (14)
[[Phi].sub.m] = 1 - x/1 + [Epsilon]x; A(T) = [C.sub.1] - [C.sub.2][(T - [T.sub.gp]).sup.2] (15)
[Mathematical Expression Omitted] (16)
where [k.sub.t0], [k.sub.p0] are the termination and propagation rates in absence of gel and glass effect, [Mathematical Expression Omitted], [Mathematical Expression Omitted], [E.sub.[Theta]p], [E.sub.[Theta]t], [C.sub.1], [C.sub.2], B are the gel and glass effect parameters, and [T.sub.gp] is the glass transition temperature.
Molecular Weight Distribution
In the classical approach to modeling free radical polymerization, the equations that describe the variation with time of the first few moments of the molecular weight distributions of the growing and of the dead polymers are commonly used in order to get average properties of the distributions and to avoid the tedious solution of a large number of equations for the concentrations of the growing and dead chains of a specified length.
The number ([M.sub.n]) and weight ([M.sub.w]) average chain lengths can be obtained from the moments as follows: [M.sub.n] = ([[Lambda].sub.1] + [[Mu].sub.1])/([[Lambda].sub.0] + [[Mu].sub.0]); [M.sub.w] = ([[Lambda].sub.2] + [[Mu].sub.2])/([[Lambda].sub.1] + [[Mu].sub.1]) and analogous relations can be derived for the variances of the distributions and other statistical properties.
This approach is particularly appropriate for continuous processes where uniform conditions are kept in the reactor over time and space. In this case, the shape of the distribution is known, and few statistical properties characterize completely the distribution itself.
Such is not the case for the process under examination. Let us consider the individual droplet first: through the course of reaction there are relevant changes in its temperature that cause the propagation and termination rates constants to change and the instantaneous molecular weight distribution (MWD) to change with them.
The change of temperature is not the only cause of variation of the kinetic rates; the gel and glass effects, which take place as reaction proceeds, further modify the instantaneous MWD so that the local cumulative MWD may be of very different shapes.
Moreover, at a fixed time, different MWD are being produced at different spatial positions along the reactor, so that global values of the instantaneous moments can be obtained only through integration and averaging over the global volume.
The first moment of the instantaneous global MWD can be computed as:
[Mathematical Expression Omitted] (17)
which for spherical symmetry becomes:
[Mathematical Expression Omitted] (18)
Similar relations can be written for the other moments as well. It is important to point out, however, that a description of a distribution curve through its moments is relevant only when the curve has a known unimodal shape, but it becomes almost insignificant for less regular curves (bimodal, trimodal, or with relevant tails). However, this happens in this case for the above-mentioned reasons (temperature variations, gel effect, glass effect, radial variations). This is the reason why it is appropriate to consider the entire distribution curves. However, instead of solving the equations in the concentrations of the growing and dead chains of all lengths of interest, integration over time and space of the local instantaneous distributions is accomplished.
In free radical polymerization, the local instantaneous number chain length distribution [f.sub.n] and weight chain length distribution [f.sub.w] of the growing radicals have the following expressions:
[f.sub.n](n, r, t) = 1/v[e.sup.-n/v]; [f.sub.w](n, r, t) = n[f.sub.n]/[integral of] n[f.sub.n]dn [between limits] [infinity] and n=1 = n/[v.sup.2] [e.sup.-n/v] (19)
where v = [[Lambda].sub.1]/[[Lambda].sub.0] is the number average chain length of the radicals. These relationships are valid if the propagation and termination rate constant do not depend on chain length and the quasi-steady-state approximation is valid for the growing radicals. It is expected, therefore, that the errors in the distributions tend to increase with the average life of the radicals (gel effect and glass effect). Under such approximations, the number average chain length of the instantaneous distribution can be computed as the ratio of the first two moments of the growing radicals. Moreover, if termination by disproportionation dominates, as in the present case, the instantaneous distributions of the dead polymer are simply equal to the radical distributions. It is apparent, however, that these distributions change over time and space as the result of the changes in [k.sub.p] and [k.sub.t].
The local cumulative distributions [F.sub.n] and [F.sub.w] can be obtained by integration over conversion and normalization of the local instantaneous distributions as follows:
[Mathematical Expression Omitted] (20)
[Mathematical Expression Omitted] (21)
It is relevant that the number chain length distribution needs to be weighted by the inverse of the number average chain length, while no weight is necessary for the weight chain length distribution.
The global cumulative distributions [Mathematical Expression Omitted] and [Mathematical Expression Omitted] can be obtained by further integration over volume and normalization:
[Mathematical Expression Omitted] (22)
[Mathematical Expression Omitted] (23)
It is important to point out once more that only the first two moments of the growing radical distributions as a function of time and space are necessary in order to obtain these distributions.
If termination by combination or chain transfer is not negligible, slightly more intricate expressions of the local instantaneous distributions may be derived, but analogous integrations to obtain the global distributions can be accomplished (36).
PARAMETRIC SENSITIVITIES: DEFINITION AND COMPUTATION
The model that results from the physical description of the system is composed of coupled partial differential equations: the material balances, which do not contain any transport term, and the thermal balance, which contains the transport term as well. The system states are therefore the concentrations, the temperature, and the first moments of the distributions as a function of both time and space. How the states of the process are sensitive to variation in the system parameters is quantitatively evaluated by introducing the parametric sensitivities.
Let [Theta] be the vector of the model parameters; the sensitivity of variable [u.sub.j](r, t) with respect to the parameter [[Theta].sub.i] is given by:
[Mathematical Expression Omitted] (24)
The computation of all parametric sensitivities involves the so-called sensitivity equations, which are derived by partial differentiation of the model equations with respect to the parameter vector [Theta].
Solution of both the system equations and the sensitivity equations is accomplished by using a software developed by Caracotsios and Stewart that extends an implicit integrator for ordinary differential-algebraic systems to deal with partial differential equations. This is obtained by finite-difference approximations of the r-derivatives on a specified spatial grid that transform the system into an ordinary differential-algebraic one. The sensitivity equations obtained by partial differentiation of the transformed system are at this stage ordinary differential ones and provide the parametric sensitivities at the grid points.
The implicit integrator DDASSL (37) employs a variable-order, variable-step predictor-corrector approach initiated by Gear (38) and is designed to handle stiff systems of coupled ordinary-differential and algebraic equations. The integrator automatically adjusts the approximation order and the step size, approximates the time derivatives by backward differences, and solves the resulting algebraic corrector equations by a modified Newton method. A direct algorithm computes at the end of each time step the parametric sensitivities.
The validation of the model has been the subject of a previous work (3). Measured and computed temperature profiles have been compared in a number of cases. As an example, Figs. 2 and 3 show the experimental values of the temperature in the positions of the four thermocouples and the results of the model for the following operating conditions: [I.sub.0] = 0.02 mol/L, R = 5 cm, [r.sub.f] = 5 (water to monomer ratio in the original feed) [T.sub.j] = 70 [degrees] C. Good qualitative agreement can be detected. However, some errors are present on the temperature profiles at short times. This is related to one of the modeling assumptions that have been chosen. The spherical symmetry of all the reactor variables implies that no heat dispersions are present on the top of the reactor. Therefore dispersions both through the upper surface and through the metallic agitator (which is not working but is still there during reaction) are neglected. Given that the conductivity parameters [[Alpha].sub.1] and [[Alpha].sub.2] have been chosen to fit the experimental data with particular care to the position in time and space where the maximum in temperature occurs, an overestimation of the gel conductivity due to the neglected dispersions take place. This overestimation explains the discrepancies in the temperature profiles at short times. However, the possible increase in model prediction ability that would take place if a model with three spatial coordinates and time as independent variables were built does not justify the enormous increase in mathematical and computational complexity. In any case, the estimation of the value of the maximum temperature and the time at which it occurs is very good for many different operating conditions (see Table 2 as described below).
Figure 4 shows the computed weight molecular weight distribution for different global conversions for the same operating conditions.
RESULTS AND DISCUSSION
The sensitivity data that are computed by the software are so voluminous that their interpretation is not straightforward. This is due to the large number of system states, of grid points and of parameters. Therefore, the most significant data need to be chosen in order to understand the sensitivity to small variations in the operating conditions, in the physical parameters and in the kinetic constants.
A good indicator of the system state is the local temperature, which is related both to the system kinetics and to the heat transfer process; indeed, any situation that can lead to runaway is connected to rapid changes in temperature. Moreover, an upper bound is imposed on temperature since the stability of the gel is guaranteed only [less than]90 [degrees] C.
The maximum temperature is reached in the center of the reactor. This result is confirmed by all the experiments and all the simulations that have been accomplished and can be easily understood from a physical analysis of the heat fluxes in the reactor. As long as the maximum temperature is not in the center of the reactor, a heat accumulation takes place as the result of both the local reaction and of a spatial heat flux. When the temperature in the center of the reactor finally overcomes the values in the surrounding area, heat dissipation is made difficult by the large distance from the refrigerating jacket.
The developed analysis therefore suggests that the temperature in the center of the reactor is a good indicator of the stability of the reactor. The parametric sensitivity of this state is used as a measure of the process sensitivity to variations in a specified parameter.
Since results relating only to temperature sensitivity [TABULAR DATA FOR TABLE 2 OMITTED] in the center of the reactor are presented, in the following the symbol S([[Theta].sub.i]) will refer to the value:
[Mathematical Expression Omitted] (25)
The parameter vector is composed of 21 elements: [I.sub.0], [T.sub.j], [r.sub.f], U, [[Alpha].sub.1], [[Alpha].sub.2], [A.sub.d], [E.sub.d], [A.sub.p], [E.sub.p], [A.sub.t], [E.sub.t], [Mathematical Expression Omitted], [Mathematical Expression Omitted], [E.sub.[Theta]p], [E.sub.[Theta]t], [C.sub.1], [C.sub.2], B, f, [Delta][H.sub.r]; the meanings of the symbols has been made clear in the development of the model.
For clarity, these parameters are divided into three groups: 1) operating conditions, 2) kinetic parameters, and 3) heat transfer parameters.
Before starting the polymerization, in order to have appropriate reaction times and desired product properties, the main parameters that must be chosen are the initiator concentration [I.sub.0], the jacket temperature [T.sub.j], and the water to monomer ratio in the initial mixture [r.sub.f].
[I.sub.0] strongly affect both the kinetics and the molecular weight of the produced polymer. [T.sub.j], which is kept constant throughout the course of reaction, has an impact on the reaction rate and therefore on internal temperature profiles; [r.sub.f] determines the heat per unit volume that is developed and affects both the thermal conductivity and the thermal capacity of the reacting mixture.
It is therefore important to understand if there is any risk in varying these parameters (e.g., to change the product properties or in the scale-up of the reactor). In this study, the base case is characterized by the following operating conditions: R = 5 cm, [I.sub.0] = 0.05 mol/L, [T.sub.j] = 70 [degrees] C; the values of all other parameters are indicated in Table 1.
Figure 5 shows both the temperature profile in the center of the reactor and the value of the parametric sensitivity with respect to the initiator concentration S([I.sub.0]). For small times, the value of S([I.sub.0]) is very low: temperature gradually increases to the jacket level, decomposition of BPO is very slow, conversion is very low, and the increase in temperature is mainly due to the heat transfer and only marginally to the heat of reaction. Above 70 [degrees] C, the reaction rate rises quickly, and, as soon as the gel effect takes place, the increase in temperature is very fast and the direction of the heat flux changes. This happens faster at higher initiator concentrations, so that the sensitivity increases significantly and reaches a high peak in a small amount of time.
Beyond the peak, there is a very fast decrease, sensitivity changes sign around the value of the temperature peak and then goes to zero for long times. This is what is expected since at long times, the overall temperature goes to the jacket value no matter how much initiator there was at the beginning.
It is not as intuitively evident that the value of the sensitivity is close to zero at the temperature peak. It is apparent from the shape of the sensitivity function that as the quantity of initiator increases, the temperature peak moves to shorter times. Given that the derivative of the temperature in the center of the reactor with respect to time is zero at the peak, such a temperature as a function of time and initiator concentration has both derivatives at the peak value close to zero.
To better understand what happens, values of the sensitivities has been computed for different values of [I.sub.0] (0.1 mol/L and 0.02 mol/L in [ILLUSTRATION FOR FIGURE 6 OMITTED]). The same qualitative behavior is observed and again a value close to zero is obtained in proximity of the temperature peak. This corresponds in the time-initiator plane to a wave of temperature moving to smaller times for high initiator concentrations but without significantly changing its maximum value.
This somewhat striking result has been confirmed by experimental results: Figs. 7 and 8 show the experimental and computed temperature, respectively, in the center of the reactor for the three different initiator concentrations. Despite the large range of variations in initiator concentration, no significant effects on the maximum temperature are detected; therefore there is no risk of overcoming the maximum temperature allowed for the stability of the gel. This interesting result can be qualitatively interpreted by considering the thermal effects related to the heat of reaction and to the dissipation through the gel. Since only heat transfer through conduction is present, it is clearly not possible to operate at isothermal conditions. However, given that the temperature rises gradually from the outside to the inside of the reactor, the reaction starts at different times and the maximum local production of heat (corresponding to the gel effect region) does not take place at the same time in the whole reactor, as would be the case for a stirred one. A decrease (increase) in the initiator concentration slows down (hastens) the moving of the temperature wave from the outside to the inside, leading to longer (shorter) reaction times. However, it does not significantly change the temperature peak, which is mainly related to the availability of monomer, to the volume involved, and to the distance from the external jacket.
In the following, an initiator concentration [I.sub.0] = 0.1 mol/L will be considered.
Interpretation of the behavior of the sensitivity to the jacket temperature S([T.sub.j]) is straightforward [ILLUSTRATION FOR FIGURE 9 OMITTED]. Higher jacket temperatures correspond to higher temperatures in the reactor, a faster reaction rate, smaller reaction times, and a higher peak. Sensitivity moves from zero as soon as the temperature wave reaches the center of the reactor, has a high peak during the gel effect, is close to one around the maximum temperature, and, after a small region with negative values, goes to one for long periods. It indicates therefore that an increase of the jacket temperature corresponds to a comparable increase of the maximum temperature; the jacket temperature must therefore be carefully chosen in order not to overcome the maximum permissible temperature. The sensitivity to the water-to-monomer weight ratio [r.sub.f] is shown in Fig. 10. As anticipated above, [r.sub.f] has two main effects on the thermal behavior of the system. First, it influences the heat transfer, as it makes heat conductivity and heat capacity vary (see the Heat Transfer Parameters section below for comparison). Second, it changes the heat that is generated per unit volume of the reactor. The two effects explain the variation from zero of S([r.sub.f]) even for short times [as it is the case for S([T.sub.j])] and its negative value throughout the course of reaction. Variations in [r.sub.f], therefore, correspond to changes not only in the reaction time but also in the reaction temperature. The value of [r.sub.f] must therefore be carefully chosen at the design stage; however, once the desired value has been determined, no uncertainties are expected at the operating stage.
Several experiments have been undertaken with variations in the operating conditions. In Table 2, the maximum temperature ([T.sub.max]) and the time ([t.sub.max]) at which it occurs obtained by experiments and by simulations are compared; also, the operating conditions, which are varied one at a time, are reported in the Table. It is important to note that the peaks in temperature are so sharp that errors of [approximately]1 [degrees] C in the estimate of the maximum value are surely acceptable. Moreover, small errors in the position of the thermocouple in the gel may give relevant variations in the maximum temperature. A decrease in the initial initiator concentration has a relevant effect on the time at which the maximum temperature occurs without significantly changing its value. Variations in the jacket temperature and in the water-to-monomer weight ratio have more relevant effects on the maximum value. In particular, when [r.sub.f] = 3, a maximum of 94 [degrees] C is obtained by the model, which corresponds to dissolution of the gel, which takes place experimentally (no experimental data are significant after the gel dissolves). All the computed and experimental results reported in Table 2 agree with the sensitivity analysis reported in this section. The last two lines of the Table refer to a variation in the agarose concentration and will be explained in the sequel.
The kinetic parameters include most of the model parameters and in particular those on which the largest uncertainty is present. For most of them, however, considerations similar to the one developed for the initiator concentration apply.
In Figs. 11 through 14, sensitivities to f, [A.sub.p] (frequency factor of the propagation kinetic constant), [E.sub.p] (activation energy of the propagation kinetic constant), [Mathematical Expression Omitted], [Mathematical Expression Omitted], [E.sub.[Theta]t], [E.sub.[Epsilon]p] and [Delta][H.sub.r] are reported. Here again, sensitivities get close to zero at the maximum temperature, so that even a strong variation in these parameters corresponds to a shift in time of the temperature peak but not to a significant variation in the peak value.
A different behavior is shown by the last two parameters. As far as [Delta][H.sub.r] is concerned, sensitivity has a maximum value in the gel effect region but does not become negative even for long times. This can be expected since a change in the heat of reaction causes a change in the same direction of the total heat produced. However, the sensitivity value is not particularly high; moreover, the heat of reaction is one of the few kinetic parameters that can be determined very precisely experimentally.
Also, the other parameter ([E.sub.[Theta]p]) with a peculiar behavior of the sensitivity is not particularly relevant; [E.sub.[Theta]p] is the activation energy for the monomer diffusion during the glass effect. A higher value of [E.sub.[Theta]p] corresponds to a higher [[Theta].sub.p] and therefore to a smaller [k.sub.p]. The reaction rate therefore gets smaller, and this corresponds to the negative value in the gel effect region. Despite the peak in this region, which is the most sensitive to all parameters, this parameter is particularly relevant when the glass effect takes place, that is to say, in the region where reaction is slowed down and no significant thermal effects are present. This explains the very small absolute values of the sensitivities to [E.sub.[Theta]p]. No significant changes in the temperature maximum are therefore expected, even for relevant errors in its estimated value.
Heat Transfer Parameters
In Figs. 15 and 16, sensitivities to the heat transfer coefficients and to the coefficient appearing in the heat conductivity expression are reported. The shape of the curves is peculiar. However, it can be explained considering that the values of these sensitivities depend mainly on two factors: the local temperature gradient (which is closely related to the difference between the local temperature and the jacket temperature) and the speed of variation of temperature. The larger these two values, the larger is the sensitivity to the heat transfer parameters. At short times, the temperature does not move from its original value so that the sensitivity is almost zero. The sensitivity gets larger than zero as soon as the temperature starts rising, and it gradually increases, reaching a maximum when it is still quite far from the jacket value and has a very high speed of variation. When the temperature gets close to the jacket value so that the conductive heat flux is quite small, the sensitivity gets smaller, since variations. in the temperature are small no matter how large the heat transfer parameters are. After that, however, the rate of reaction increases at a fast rate so that both temperature and its speed of variation increase rapidly and the peak for the sensitivity is reached. A steep drop is present, as for the other parameters around the temperature maximum, where the sensitivity value is close to zero.
Also in this case, therefore, variations in the values of U, [[Alpha].sub.1], [[Alpha].sub.2] (or errors in their estimate) can lead to errors in the expected reaction time but not to relevant ones in the estimate of the maximum temperature achieved in the reactor. This is due to the behavior of the sensitivity in the time interval around the maximum temperature; indeed, this behavior resembles the one of the sensitivities to most of the other parameters that have been analyzed so far. The variations both in the maximum temperature and in the time at which this maximum occurs have been tested with experiments in which the agarose concentration (and therefore the gel conductivity) has been varied. A decrease in the agarose concentration corresponds to a less compact gel and therefore to a higher conductivity; the maximum temperature decreases but not significantly (Table 2), and the reaction takes place in a shorter time, as expected. It is important to note that when the agarose concentration is varied, the conductivity parameters [[Alpha].sub.1] and [[Alpha].sub.2] are recomputed in order to fit the experimental data (% agarose = 0.67, [[Alpha].sub.1] = 2, [[Alpha].sub.2] = 6.6 x [10.sup.-4] [K.sup.-2]; % agarose = 0.5, [[Alpha].sub.1] = 2, [[Alpha].sub.2] = 7.2 x [10.sup.-4] [K.sup.-2]). A good agreement has been obtained and the predictions resulting from the sensitivity analysis have been confirmed.
A detailed analysis for a new suspension polymerization process has been presented. The process is characterized by peculiar temperature profiles given that the polymerizing particles are embedded into a gel phase that prevents heat dissipation through convection; this is done so that no suspending agent is needed to avoid coalescence, therefore, high purity polymers are produced.
A first principle model has been used to study MMA polymerization in a gelled phase of water and agarose. It provides a good description of the temperature profiles that have been obtained through experiments in a lab reactor.
Given that different thermal histories take place at different radial positions and that a constraint on the maximum temperature that guarantees the stability of the gel is present, parametric sensitivity of the maximum temperature in the reactor is particularly relevant. Theoretical and experimental sensitivity studies give necessary information for reactor scale-up and process modifications and also help understanding the physics and the kinetics of the process. Moreover, they can be used to predict any risk of thermal runaway.
The operating conditions have first been accounted for; changes in initiator concentration lead to a relevant change in the time at which the maximum temperature is achieved, but only to a slight variation in the maximum. Significant variations in the maximum temperature are related to variations in the jacket temperature. The effects of changes in the water to MMA weight ratio in the initial mixture are even more relevant; this parameter therefore needs to be carefully chosen. These results have been obtained theoretically by means of the mathematical model; moreover, they have been confirmed by several experiments.
Subsequently, the analysis has been extended to the kinetic and heat transfer parameters for which only the theoretical study is meaningful. Almost all kinetic and heat transfer parameters have an effect, which is qualitatively similar to the initiator concentration. The few parameters that behave differently either do not have relevant uncertainties or do not have a large effect on the temperature peak. Therefore, no risk of runaway or of instability of the gel has been detected once the jacket temperature and the water to MMA weight ratio have been appropriately chosen.
The developed process is therefore interesting for wider applications and can be extended to the polymerization of other monomers for which high purity products are desired.
Giovanni Polacco thanks Scuola Normale Superiore and Montecatini SpA for supporting his participation in this work by a PhD grant.
A = Frequency factor.
B, [C.sub.1], [C.sub.2] = Gel and glass effect parameters.
[c.sub.p] = Specific heat.
[Delta][H.sub.r] = Heat of reaction.
E = Activation energy.
[E.sub.[Theta]p] = Glass effect parameter.
[E.sub.[Theta]t] = Gel effect parameter.
f = Initiator efficiency.
[f.sub.n](n) = Local instantaneous number chain length distribution.
[f.sub.w](n) = Local instantaneous weight chain length distribution.
[F.sub.n](n) = Local cumulative number chain length distribution.
[F.sub.w](n) = Local cumulative weight chain length distribution.
[Mathematical Expression Omitted] = Global cumulative number chain length distribution.
[Mathematical Expression Omitted] = Global cumulative weight chain length distribution.
I = Initiator concentration.
K = Thermal conductivity.
[k.sub.d] = Rate of initiator decomposition.
[k.sub.i] = Rate of initiation.
[k.sub.p] = Rate of propagation.
[k.sub.p0] = Rate of propagation in absence of gel and glass effect.
[k.sub.t] = Rate of termination.
[k.sub.to] = Rate of termination in absence of gel and glass effect.
[k.sub.tc] = Rate of termination by combination.
[k.sub.td] = Rate of termination by disproportionation.
M = Monomer concentration.
[M.sub.n] = Number average chain length.
[M.sub.w] = Weight average chain length.
[P[prime].sub.n] = Growing radical of length n.
[P.sub.n] = Dead polymer of length n.
R[prime] = Primary radical.
R = Reactor radius.
[r.sub.f] = Water to monomer ratio in the initial mixture.
S = Sensitivity.
T = Temperature.
t = Time.
[T.sub.gp] = Glass transition temperature.
U = Heat transfer coefficient for the reactor jacket.
u = System state vector.
V = Reactor volume.
x = Monomer conversion.
[[Alpha].sub.1], [[Alpha].sub.2] = Parameters for the thermal conductivity.
[Epsilon] = Volume expansion factor.
[[Epsilon].sub.[Sigma]] = Volume fraction of the reacting phase.
[Mathematical Expression Omitted] = Glass effect parameter.
[Mathematical Expression Omitted] = Gel effect parameter.
[[Lambda].sub.i] = i-order moment of the growing chains.
[[Mu].sub.i] = i-order moment of the dead polymer.
[Rho] = Density.
[Theta] = Parameter vector.
0 = Initial conditions.
c = Computed value.
d = Initiator decomposition.
e = Experimental value.
m = Monomer.
p = Polymer/propagation.
t = Termination.
1. P. Giusti, M. Mantilli, M. Palla, and G. Pizzirani, Italian Patent 67967 A90 (1990).
2. P. Giusti, L. Lazzeri, G. Pizzirani, G. Polacco, C. Rizzo, and M. Palla, J. Mater. Sci.: Mater. Med., 5, 858 (1994).
3. G. Polacco, D. Semino, and C. Rizzo, J. Mater. Sci.: Mater. Med., 5, 587 (1994).
4. M. Morbidelli and A. Varma, Chem. Eng. Sci., 40, 2165 (1985).
5. H. Rabitz, M. Kramer, and D. Dacol, Ann. Rev. Phys. Chem., 34, 419 (1983).
6. J. W. Tilden, V. Constanza, G. J. McRae, and J. H. Seinfeld, in Modeling of Chemical Reaction Systems, p. 69, K. H. Ebert, P. Denfhlart, and W. Jager, eds. SpringerVerlag, Berlin (1981).
7. N. N. Semenov, Z. Phys., 48, 571 (1928).
8. J. B. Rawlings and W. H. Ray, Polym. Eng. Sci., 28, 237 (1988).
9. M. Morbidelli and A. Varma, Chem. Eng. Sci., 43, 91 (1986).
10. B. Kapoor, S. K. Gupta, and A. Varma, Polym. Eng. Sci., 29, 1246 (1989).
11. M. Tjahjadi, S. K. Gupta, M. Morbidelli, and A. Varma, Chem. Eng. Sci., 42, 2385 (1987).
12. M. Caracotsios and W. E. Stewart, Computers Chem. Eng., 9, 359 (1985).
13. S. Zhu and A. E. Hamielec, Polymer, 32, 3021 (1991).
14. S. Zhu and A. E. Hamielec, Macromolecules, 22, 3098 (1989).
15. W. Y. Chiu, G. M. Carratt, and D. S. Soong, Macromolecules, 16, 348 (1983).
16. P. E. Baillagou and D. S. Soong, Polym. Eng. Sci., 25, 213 (1985).
17. E. H. Trommsdorf, H. Koenle, and P. Lagally, Makromol. Chem., 1, 169 (1948).
18. R. G. W. Norrish and E. F. Brookman, Proc. Roy. Soc. (London), A171, 147 (1939).
19. G. W. Shultz and G. Harbort, Makromol. Chem., 1, 106 (1947).
20. A. W. Hui and A. E. Hamielec, J. Appl. Polym. Sci., 16, 749 (1972).
21. N. Friis and A. E. Hamielec, J. Polym. Sci., (A-1), 12, 251 (1974).
22. R. T. Ross and R. L. Laurence, AIChE Symp. Ser., 160, 74 (1976).
23. A. Husain and A. E. Hamielec, J. Appl. Polym. Sci., 22, 1207 (1978).
24. C. C. Lin and Y. F. Wang, J. Appl. Polym. Sci., 26, 3909 (1981).
25. J. N. Cardenas and K. F. O'Driscoll, J. Polym. Sci., Polym. Chem. Ed., 14, 883 (1976).
26. S. K. Soh and D. C. Sundberg, J. Polym. Sci., Polym. Chem. Ed., 20, 1315 (1982).
27. A. D. Schmidt and W. H. Ray, Chem. Eng. Sci., 36, 1401 (1981).
28. T. J. Tulig and M. V. Tirrel, Macromolecules, 14, 1501 (1981).
29. K. F. O'Driscoll, Pure Appl. Chem., 53, 617 (1981).
30. A. E. Hamielec, Chem. Eng. Commun., 24, 1 (1983).
31. V. R. Kumar and S. K. Gupta, Polymer, 32, 323 (1991).
32. N. R. Vaid and S. K. Gupta, Polym. Eng. Sci., 31, 1708 (1991).
33. J. S. Chang and W. J. Chen, Canadian J. Chem. Eng., 70, 523 (1992).
34. P. E. Baillagou and D. S. Soong, Chem. Eng. Sci., 40, 75 (1984).
35. D. Achilias and C. Kiparissides, J. Appl. Polym. Sci., 35, 1303 (1988).
36. G. Polacco and D. Semino, submitted to Makromol Chem., Theory and Simulations.
37. L. R. Petzold, in Scientific Computing, R. S. Stepleman et al., eds., North-Holland, Amsterdam (1983).
38. C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, N.J. (1971).
|Printer friendly Cite/link Email Feedback|
|Author:||Polacco, Giovanni; Semino, Daniele; Palla, Maurizio|
|Publication:||Polymer Engineering and Science|
|Date:||Aug 1, 1996|
|Previous Article:||The glass transition temperature of nitrated polystyrene/poly(acrylic acid) blends.|
|Next Article:||The relative influences of molecular structure on brittle fracture by fatigue and under constant load in polyethylenes.|