# Simulation of bulk free radical polymerization of styrene in tubular reactors.

INTRODUCTIONPolystyrene is usually produced by free radical bulk polymerization of styrene in a batch tank reactor (expandable polystyrene) or continuous stirred tank reactor (general purpose and high impact polystyrene). The use of tubular reactors is being advocated by many researchers due to their simplicity and low cost. However, because of the need of the quick removal of the exothermic heat of reaction and the possibility of tubes plugging due to the high viscosity of the product, static mixers and loop reactors have been suggested to improve the performance of tubular reactors (1).

Previous investigation included some experimental work, Wallis, et al. (2) and Valsamis and Biesenberger (3) and numerous theoretical works on modeling and simulation of polystyrene tubular reactors (2, 4-8). Because of the large dependence of the viscosity of the reaction mixture on temperature, composition and molecular weight, there is a strong interaction between the fluid mechanics modeling equations and the kinetics modeling equations. Early work (2, 5) neglected this interaction by assuming a parabolic profile for the velocity. Later work, (4, 6, 8) included this interaction with different degrees of rigor. The most rigorous model used was that of Stevens and Ray (7). In this paper, however, we use the model presented by Chen and Nauman (8) due to its simplicity and ability to represent the interaction by taking into consideration radial and axial velocity change due to viscosity change which in turn depends on temperature and conversion.

A stability consideration has limited the tube radius to 2 cm to avoid velocity profile elongation and thermal runaway (6).

Sensitivity studies (2, 5) indicated that an increase in inlet temperature, tube wall temperature, and tube radius will lead to a decrease in average molecular weight and an increase in polydispersity; and thus to a poor quality product.

The assumption of a parabolic velocity profile is an over simplification since the tube wall viscosity is usually much greater than the tube center viscosity. Taking radial viscosity changes into consideration to calculate axial velocity leads to a system of equations that are difficult to solve. Vrentas and Huang (9) found that a fine finite difference mesh in the axial and radial directions is needed to obtain accurate solution for high viscosity changes. As many as 480 radial increments and 2000 axial increments are required in extreme cases. Agarwal and Kleinstreuer (6) used 20 and 100 increments in the radial and axial directions respectively. The use of adaptive grids will therefore be recommended to solve the system equations. Stevens and Ray (7) used B-splines on adaptive grids to solve their rigorous model.

Alternatively the method of lines can be used by discretizing the equation in the radial direction using finite difference, finite elements or collocation methods. The system equations become ordinary differential equations of the initial value type. The system equations are usually stiff and the Gear method is thus recommended for their solution. Chen and Nauman (8) used ten finite difference mesh in the radial direction.

The purpose of this work is to investigate the numerical solution of the modeling equations using the method of orthogonal collocation and to suggest design and operation modifications to improve the performance of tubular reactors. Specifically we investigate the spatial changes in wall temperature and tube radius.

MODEL

Assumptions used for simplification of model are:

1. The flow is laminar, steady, and axisymmetric.

2. Viscosity depends on temperature and composition but independent of shear rate.

3. The velocity is fully developed in the axial direction.

4. Density, viscosity, velocities, component concentration, temperature, and molecular weights vary both radially and axially. Other properties are constant.

5. The monomer obeys Fick's law of diffusion. Polymer is non-diffusing, but overall continuity is satisfied by equal mass counter-displacement.

6. Mass diffusion and heat conduction in the axial direction are negligible.

7. Viscous heating effects are negligible.

8. The reactor wall temperature is constant.

These assumptions were used by Chen and Nauman (8). Transport equations in cylindrical coordinates based on above assumptions are:

1. Overall mass balance

[Mathematical Expression Omitted]

or in differential form

1/r [Delta]/[Delta]r ([Rho]r[v.sub.r]) + [Delta]/[Delta]z ([Rho][v.sub.z]) = 0 (2)

2. Momentum balance

The z-component equation of motion is

0 = -[Delta]P/[Delta]z + 1/r [Delta]/[Delta]r ([Mu]r [Delta][v.sub.z]/[Delta]r) (3)

The boundary conditions are

[Delta][v.sub.z]/[Delta]r= 0 at r = o (4a)

[v.sub.z] = 0 at r = R (4b)

By integration twice

[v.sub.z](r, z) = 1/2 (-dP/dz) [integral of] r/[Mu](r, z) dr between limits R and r (5)

Substitute in Eq 1 the value of [v.sub.z](r, z)

(-dP/dz) = [Q.sub.m]/[Pi] 1/[integral of] r[[Rho].sub.r, z] between limits R and o [integral of] r/[Mu](r, z) dr dr between limits R and r (6)

and also

[v.sub.z](r, z) = [Q.sub.m]/2[Pi] [integral of] r/[Mu](r, z) dr between limits R and r/[integral of] r[Rho](r, z) between limits R and o [integral of] r/[Mu](r, z) dr dr between limits R and r (7)

The radial velocity is then

[v.sub.r](r, z) = -1/r[Rho](r, z) [Delta]/[Delta]z [integral of] r[Rho](r, z) [v.sub.z](r, z) dr between limits r and o (8)

3. Component Mass Balance

The differential mass balance equation is given by

[Delta]/[Delta]z ([Rho][v.sub.z][w.sub.m]) + 1/r [Delta]/[Delta]r ([Rho]r[v.sub.r][w.sub.m]) = [D.sub.m]/r [Delta]/[Delta]r ([Rho]r [Delta][w.sub.m]/[Delta]r) - [R.sub.p] (9)

The boundary conditions are

[w.sub.m] = [w.sub.m, in] at z = 0 (10a)

[Delta][w.sub.m]/[Delta]r = 0 at r = 0 (10b)

[Delta][w.sub.m]/[Delta]r = 0 at r = R (10c)

4. Energy Balance

The differential energy balance equation is given by:

[C.sub.p] [Delta]/[Delta]z ([Rho][v.sub.z]T) + [C.sub.p] 1/r [Delta]/[Delta]r ([Rho]r[v.sub.r]T) = K([[Delta].sup.2]T/[Delta][r.sup.2] + 1/r [Delta]T/[Delta]r) + (-[Delta]H) [R.sub.p] (11)

The boundary conditions are

T = [T.sub.in]. at z = 0 (12a)

[Delta]T/[Delta]r = 0 at r = 0 (12b)

T = [T.sub.w] at r = R (12c)

Kinetics and Physical Properties

The detailed reaction mechanism is given by Chen and Nauman (8) for the thermally initiated free radical polymerization of styrene. The kinetic parameters and physical property data used are given in Table 1.

Molecular Weight Determination

The instantaneous average molecular weights are given by

[Mathematical Expression Omitted]

[Mathematical Expression Omitted]

Table 1. Kinetic Parameters and Physical Property Data.

[K.sub.i] = 2.019 x [10.sup.1] exp(-13,810/T), [m.sup.6]/([kg.sup.2])(s) [K.sub.p] = 1.009 x [10.sup.5] exp(-3557/T), [m.sup.3]/(kg)(s) [([K.sub.tc]).sub.0] = 1.205 x [10.sup.7] exp(-844/T), [m.sup.3]/(kg)(s) [Mathematical Expression Omitted] [A.sub.1] = 2.57 - 5.05 x [10.sup.-3]T [A.sub.2] = 9.56 - 1.76 x [10.sup.-2]T [A.sub.3] = 3.03 + 7.85 x [10.sup.-3]T [K.sub.tr, m] = 2.218 x [10.sup.4] exp(-6377/T), [m.sup.3]/(kg)(s) [D.sub.m] = 2 x [10.sup.-9] [m.sup.2]/s [C.sub.p] = 1.88 x [10.sup.3] J/(kg)([degrees] K) K = 0.126 J/(m)(s)([degrees] K) [Delta]H = -6.70 x [10.sup.5] J/kg [Rho] = (1174.7 - 0.918T)(1 - [w.sub.p]) + (1250.0 - 0.605T) [w.sub.p], kg/[m.sup.3] [Mu] = exp{ - 13.04 + 2013/T + M[w.sup.0.18] [Mathematical Expression Omitted]

where

[[Mu].sub.0] = [R.sub.p] ([C.sub.m] + [Beta]/2)

[[Mu].sub.1] = [R.sub.p]([C.sub.m] + [Beta] + 1) [is approximately equal to] Rp

[[Mu].sub.2] = [R.sub.p](2[C.sub.m] + 3[Beta])/[([C.sub.m] + [Beta]).sup.2], [Mathematical Expression Omitted], [C.sub.m] = [K.sub.tr, m]/[K.sub.p]

The observed monomer chain transfer is written as

([C.sub.m]) apparent = ([C.sub.m]) actual + [B.sub.i][w.sub.p]

[B.sub.1] = -9.85 x [10.sup.-2] + 7.77 x [10.sup.-4]T - 2.04 x [10.sup.-6] [T.sup.2] + 1.79 x [10.sup.-9] [T.sup.3] (15)

The stream line average molecular weights for I axial increments are given by

[Mathematical Expression Omitted]

[Mathematical Expression Omitted]

where [w.sub.i] is the weight of polymer of molecular weight [M.sub.i]. The variable, [[Rho].sub.i + z] + [w.sub.p, i + z], represents the weight of polymers produced at axial position i + 1 after monomer reaction without diffusion in the axial increment i to i + 1. The value of ([[Rho].sub.i][w.sub.p, i + 1] - [[Rho].sub.i][w.sub.p, i]) represents the amount of dead polymer formed by reaction in the axial increment i to i + 1.

The cup-averaged molecular weights for polymer leaving the reactor is

[Mathematical Expression Omitted]

[Mathematical Expression Omitted]

METHOD OF SOLUTION

The method of orthogonal collocation (10) is used to discretize the modeling equations in the radial direction. The equations become ordinary differential equations in the axial direction which are solved by the DGEAR subroutine of the IMSL library which is based on the Gear method for the solution of stiff differential equations. In fact as the number of collocation points increases, the stiffness of the equations increases and the equations become more difficult to solve and in some cases the method fails to converge to the solution. For this purpose, we would like to discretize the equations in the radial direction with a minimum of collocation points. Convergence for the results to the solution is ascertained by their constancy as the number of collocation points increases.

One reason for the need of numerous collocation points is the steepness of the velocity profile. Thus we investigate methods which could calculate the velocity profile with the minimum of collocation points.

From Eq 5, we could use the trapizoidal rule or Simpson's rule to evaluate the integrals. The trapizoidal rule was more accurate in evaluating the integrals in Eq 5 but less accurate than direct application of the orthogonal collocation method to the differential equation (Eq 2) and the Radau quadrature to Eq 1. This last method is thus used for calculating the velocity profile in this work.

However the use of Eqs 1 and 2 to calculate velocity profile still requires a large number of collocation points.

The following transformation

u = [square root of] [v.sub.z]

has been tried and the differential equation is solved for u but this transformation did not reduce the number of collocation points for a certain accuracy.

The solution of the first order differential equation obtained by integrating Eq 2 once using collocation gives less accurate results.

For an inlet flow rate of 0.4842 kg/h, an inlet temperature of 413 K and a wall temperature of 408 K, six points collocation are enough to give an accurate solution for a radius of 0.003 m giving a maximum dimensionless axial velocity of 2.72. Eight points collocation gives accurate solution for a radius of 0.004 m with a maximum dimensionless axial velocity of 3.88. As the tube radius increases, velocity elongation increases and this requires more collocation points making equations stiffer and failing to converge to the solution. For the purpose of this study we limit ourselves with small radius tubes and conclude that we should be suspicious of the results obtained by using large radius tubes and only ten points of discretization in the radial direction, (4, 8). Although temperature and concentration profiles are not so sensitive to the number of collocation points, any slight change in their values causes a large change in the viscosity and hence the velocity profile.

BASE CASE

The following parameters are used for the base case study, tube radius = 0.003 m, wall temperature = 408 K, inlet temperature = 413 K, mass flowrate = 0.4842 kg/h, and tube length = 6.4 m.

Figure 1 illustrates the variation of conversion along the radial direction at different axial positions. The variations between the conversion near the wall and the tube center along the axial direction are considerable. This leads to nonuniformity of the polymer products. Due to the small conversion in this case, and lower wall temperature than inlet temperature, the temperature profile decreases from an inlet temperature of 413 K to 408.9 K at the tube center. No hot spot is observed in this case.

Figure 2 illustrates the variation of both cup-averaged weight average molecular weight (C.A.W.A.M.W.) (curve 1) and cup averaged number average molecular weight (C.A.N.A.M.W.) (curve 2), along the axial direction.

PARAMETRIC STUDY

As mentioned in the Introduction, the effect of different parameters on the product quality has been studied by previous investigators. However for comparison purposes with the results of the present study, we give in Table 2, the results of a single variation in the inlet temperature, wall temperature, inlet flowrate, tube radius for the case of constant inlet flowrate, and tube radius for the case of constant inlet mass velocity.

The leading factor affecting concentration and temperature profiles is the presence of a core in the center having a higher temperature and smaller residence time than the layer near the wall. Other factors of importance are the radial thermal conduction and the radial diffusion and radial viscosity variation with temperature and composition.

In Table 2 it is noticed that the change in the inlet temperature affects conversion very slightly because of the small tube radius used leading to the quick removal of heat of reaction. For wall temperature increase, the average reactor temperature increases leading to higher rate of polymerization, thus conversion increases, radial variation in the conversion increases, and cup averaged molecular weights decrease. The increase in the flowrate and tube radius cause also an increase in the radial variation of the conversion. A decrease in the inlet flowrate leads to higher conversion because of larger residence time. Increasing the tube radius while keeping the inlet flowrate constant causes an increase in the average conversion because less heat of reaction will be removed due to larger tube radius with higher center temperature leading to more radial variation in the conversion.

EFFECT OF WALL TEMPERATURE CHANGE

In all cases studied in this paper when one variable is changed, the others are kept at the following values of:

Tube radius = 0.003 m, wall temperature = 408 K, inlet temperature = 413 K, mass flow rate = 0.4842 kg/h, tube length = 6.4.

It is required to reduce the wall temperature in the axial direction in order to reduce the conversion near the wall resulting in reduction of molecular weight distribution. A preliminary optimization using two wall temperatures strategy has been carried out by Nauman and Mallikarjun (11). Five cases have been selected for the study; these are:

* Case 1, 2, 3; where the wall temperature decreased 3 [degrees] C every 80 cm from the tube length beginning from 403, 408, and 413 K respectively.

* Case 4, where the wall temperature starts from 387 K and increases 3 [degrees] C every 80 cm from the tube length.

* Case 5, which is identical to case 2 but the decrease is 2 [degrees] C instead.

The results of calculation are summarized in Table 3.

By comparing the first three cases, it can be noticed that by using the wall temperature variation, the outlet conversion both at the center and near the wall will decrease; meanwhile the difference between them will be reduced and this tends to reduce the molecular weight distribution. Also the cup-averaged molecular weights will increase and so will the polydispersity.

By comparing case 2 with case 4 it can be noticed that the difference between the outlet conversion at the center and the wall of the tube is higher in case 4 than in case 2, and the same applies to the cup-averaged molecular weights.

By comparing case 2 with case 5 it can be noticed that the outlet conversion both in the center and near the wall are greater in case 5 than in case 2 and the differences are greater than in the previous cases. Hence as the wall temperature is reduced in the axial direction the cup-averaged molecular weights and the polydispersity will increase.

To confirm that these conclusions are generally valid, a case of high conversion is studied. In this case the inlet temperature is increased to 463 K and the tube radius is increased to 0.0075 m. The results of the simulation are presented in the first row of Table 4.

Next the wall temperature is reduced from 408 K at the inlet by 1.5 K every 80 cm. The results of the simulation are presented in the second row of Table 4.

EFFECT OF RADIUS CHANGE (CONICAL SHAPE TUBE)

The gradual increase in the tube radius will result in velocity elongation reduction (caused by the highly increasing viscosity near the wall more than that at the center).

TABULAR DATA OMITTED

TABULAR DATA OMITTED

TABULAR DATA OMITTED

Four cases have been selected, the first two cases are tube radius increase while the third and fourth are tube radius reduction as shown in Table 5. The results are summarized in Table 6.

From cases 1 and 2 it can be noticed that with the increase in radius the outlet conversion at the center and near the wall will increase and the difference between them will decrease. Also the cup-averaged molecular weight and the polydispersity will increase.

By comparing cases 1 and 3 it can be noticed that the decrease in tube radius along the axial direction will end with large conversion at the center and near the wall, with small differences between them and greater cup-averaged molecular weight and polydispersity.

Again a high conversion case is studied. In this case the tube radius is reduced linearly along the longitudinal direction from 0.01006 to 0.0075 m. The results of the simulation are presented in the third row of Table 4.

EFFECT OF COMBINATION OF WALL TEMPERATURE CHANGE AND TUBE RADIUS CHANGE

The product properties are improved by the gradual decrease of wall temperature. Meanwhile starting with relatively large tube radius at the inlet and ending with a small radius at the outlet will result in large conversion at the center and close to that at the wall. Both effects, wall temperature change and tube radius change, will be considered simultaneously.

Based on the previous cases, two new cases are attempted;

Case 1: Case 2 in wall temperature decrease plus case 2 in radius increase.

Case 2: Case 2 in wall temperature decrease plus case 4 in radius decrease.

From comparison of the results in Table 7 with the corresponding cases in Table 6, case 2 and case 4 but without a decrease in wall temperature, it is noticed that the outlet conversion at the center and at the wall and the difference in their values are smaller in Table 7. The cup averaged molecular weights and the polydispersity are larger in Table 7. On the other hand, comparing the results in Table 7 with the corresponding case in Table 3 (case 2) with constant radius, it is noticed that the outlet conversion at the center and at the wall are larger in Table 7 but the difference in their values is smaller. The cup averaged molecular weights and the polydispersity are larger in Table 7.

Table 5. Studied Cases for Tube Radius Change. No. Tube Radius Outlet Radius 1 0.003 + 0.00003 * Z, R[R.sub.0] = 4.92 X [10.sup.-3] m 2 0.003 + 0.00004 * Z, R[R.sub.0] = 5.560 X [10.sup.-3] m 3 0.00492 - 0.00003 * Z, R[R.sub.0] = 3.00 X [10.sup.-3] m 4 0.00566 - 0.00004 * Z, R[R.sub.0] = 3.00 X [10.sup.-3] m Where Z is the axial distance in meters.

Similar observations are noticed in the last row of Table 4 giving the results for the simultaneous decrease in wall temperature by 1.5 [degrees] C/80 cm and radius change from 0.01006 m inlet diameter to 0.003 m outlet diameter for the high conversion case.

The main conclusion is that the simultaneous decrease in wall temperature and tube radius helps to reduce the difference in conversion between the center and tube wall.

CONCLUSIONS

This study indicated the usefulness of using tapered tubes with decreasing radius along the direction of the flowing reaction system as well as decreasing the wall temperature. This may be difficult to implement, but a two section reactor with different tube radii and different wall temperatures could be the practical way of implementing the findings of this work. Numerical difficulties are encountered due to the steepness of the velocity profile. Future work could be directed towards finding the optimal wall temperature and variable radius strategies. The numerical problems indicate the importance of using adaptive grids in simulating tubular polymerization reactors.

NOMENCLATURE

[A.sub.1], [A.sub.2], [A.sub.3] = Gell effect coefficients.

[B.sub.1] = Chain transfer concentration correction factor.

[C.sub.m] = Chain transfer constant for monomer.

[C.sub.p] = Heat capacity.

[D.sub.m] = Molecular diffusion coefficient for monomer.

(-[Delta]H) = Heat of reaction.

K = Thermal conductivity.

[([K.sub.tc]).sub.0] = Rate constant for combination termination at zero conversion.

M = Molecular weights.

[M.sub.n] = Cup-averaged number average molecular weight.

[M.sub.w] = Cup-averaged weight average molecular weight.

[M.sub.ns] = Number average molecular weight for streamline at the reactor exit.

[M.sub.ws] = Weight average molecular weight for stream line at the reactor exit.

[M.sub.n], i = Number average molecular weight at axial position i at any streamline.

[M.sub.w], i = Weight average molecular weight at axial position i at any streamline.

[M.sub.n], inst = Instantaneous weight average molecular weight.

[M.sub.w], inst = Instantaneous weight average molecular weight.

P = Pressure in the tube.

[Q.sub.m] = Mass flow rate.

[Q.sub.v] = Volumetric flow rate.

R = Tube radius.

[R.sub.p] = Rate of polymerization.

r = Radial coordinate.

T = Temperature.

v = Velocity.

[w.sub.p] = Polymer weight fraction.

[w.sub.m] = Monomer weight fraction.

Greek

[Rho] = Density.

[Mu] = Viscosity.

[[Mu].sub.0] = Zero polymer moment.

[[Mu].sub.1], [[Mu].sub.2] = First and second polymer moment.

[[Gamma].sub.0] = Zero radial polymer moment.

[[Gamma].sub.1], [[Gamma].sub.2] = First and second radial polymer moment.

Subscript

Z = Axial direction.

r = Radial direction.

in = Inlet conditions.

W = Wall conditions.

Superscript

- = Average value.

TABULAR DATA OMITTED

TABULAR DATA OMITTED

REFERENCES

1. N. K. Tien, F. Streiff, E. Flaschel, and A. Renken, Chem. Eng. Technol., 13, 214 (1990).

2. J. P. A. Wallis, R. A. Ritter, and H. Andere, AIChE J., 21, 686 (1975).

3. L. Valsamis and J. A. Biesenberger, AIChE Symp. Series, 72, 18 (1976).

4. C. E. Wyman and L. F. Carter, AIChE Symp. Series, 72, 1 (1976).

5. A. Husain and A. E. Hamielec, AIChE Symp, Series, 72, 111 (1976).

6. S. S. Agarwal and C. Kleinstreuer, Chem. Eng. Sci., 41, 3101 (1986).

7. C. J. Stevens and W. H. Ray, "The Mathematical Modeling of Bulk and Solution Polymerization in a Tubular Reactor," International Symposium on Computer Applications in Applied Polymer Science, 3rd Chemical Engineering Congress of North America and 195th ACS National Meeting, Toronto (1988).

8. C. C. Chen and E. B. Nauman, Chem. Eng. Sci., 44tl, 179 (1989).

9. J. S. Vrentas and W. J. Huang, Chem. Eng. Sci., 41, 2041 (1986).

10. J. Villadsen and M. L. Michelsen, Solution of Differential Equation Models by Polynomial Approximation Prentice-Hall, Englewood Cliffs, N.J. (1978).

11. E. B. Nauman and R. Mallikarjun, ACS Symp. Series 237, 305 (1984).

Printer friendly Cite/link Email Feedback | |

Author: | Soliman, M.A.; Aljarboa, T.; Alahmad, M. |
---|---|

Publication: | Polymer Engineering and Science |

Date: | Oct 15, 1994 |

Words: | 4052 |

Previous Article: | Shear and elongational behavior of linear low-density and low-density polyethylene blends from capillary rheometry. |

Next Article: | Single screw extrusion of viscoplastic fluids subject to different slip coefficients at screw and barrel surfaces. |

Topics: |