- Email: [email protected]
Dynamic optimisation of styrene polymerisation in batch reactors
European Symposiumon Computer-AidedProcess Engineering- 14 A. Barbosa-P6voa and H. Matos (Editors) 9 2004 ElsevierB.V. All rights reserved.
649
Dynamic Optimisation of Styrene Polymerisation in Batch Reactors E.E. Ekpo and I.M. Mujtaba* School of Engineering Design & Technology University of Bradford, West Yorkshire BD7 1DP, UK.
Abstract This paper focuses on the optimisation of free radical polymerisation of styrene using 2,2' azobisisobutyronitrile catalyst, (AIBN) as initiator in a batch reactor. A dynamic optimisation method is used to find the optimal temperature profile that will yield a desired level of monomer conversion and number average molecular weight in minimum time. The batch time is divided into a finite number of intervals and piecewise constant temperature is used in each interval. In each interval, the temperature and the length of the interval are optimised.
Keywords:
Batch Polymerisation, Molecular Weight Distribution, Optimal Reactor Temperature, Dynamic Optimisation
1. Introduction In today's marketplace where demands can oscillate widely in short periods, versatile systems that can adapt to changing trends are needed. This is where a batch reactor has advantage over a continuous reactor. Small batch reactors require less auxiliary equipment, and their control systems are less elaborate and costly (Ray, 1986). Polystyrene is a widely used polymer with many day-to-day uses from flimsy foam packaging to more durable plastic parts used in automobiles. It is certainly in the interest of this billion-Euro-a-year business that production be optimised. In this work, a minimum time dynamic optimisation problem is formulated for free radical bulk polymerisation in a batch reactor using temperature or coolant flow as the control variable. Here, the problem is posed as a Nonlinear Programming (NLP) problem, which is solved using an SQP-based optimisation technique. The technique has been widely used in the past (Aziz and Mujtaba, 2002). The goal is to find the optimal temperature or coolant flow profiles and initial amount of initiator (I0) to yield specified values of monomer conversion (m*) and the Number Average Chain length (Xn) in minimum time. Many chemical engineering systems are transient and are usually modelled by combinations of Differential-Algebraic Equations (DAEs) of varying complexities. Suitable numerical techniques to simulate and optimise the operation of these systems have been sought in the past. Some of these are: a. The use of Dynamic Programming techniques. b. The use of NLP based Dynamic Optimisation techniques. c. The use of Pontryagin's Maximum Principle. * all correspondences to Dr. I.M. Mujtaba, email: [email protected]
650 The last two techniques are of interest to us. Raja (1995) used the method of Lagrange multipliers to estimate the initial Initiator Concentration, I0 needed and the time it would take for the polymerisation to go on to completion and then, Pontryagin's Maximum Principle to optimise the values of X, and m (a two-step process). Here, we resort to using a Control Vector Parameterisation (CVP) technique to pose the dynamic optimisation problem as NLP problem which is then solved with SQP-based optimisation routines. The main attraction of the CVP technique is its ability to handle large systems without needing to solve an excessively large optimisation problem. CVP methods are also faster than IDP methods (Fikar et.al.,1998). As a result, the need for a two-part process as in the case of the maximum principle is removed. We also get closer values for the specified parameters in shorter batch times.
2. Model Equations The model used in this work has been developed from two sources. (Wallis et al, 1975; Young and Lovell, 1991). It is a set of DAEs represented in terms of the first three moments of the dead polymer that describe the changes taking place in the batch reactor. Monomer conversion (m) and number average molecular weight (Xn) have to be evaluated from this system of DAEs. From these, the MWD can be represented by the Shulz formula: W(j) = 10 * (-log W )b+2 jb+l Tj / r(b+2)
(1)
N(j) = (W(j)/j)/Z(W(j)/(j)
(2)
where b = -(2-PD)/(1-PD) and T = exp (-(b+ 1)/Xn). W(j) and N(j) are defined as the continuous weight and number average molecular weights respectively. Some of the model equations are presented below: d__.~! = kl (l-C) ~ (1-.~!) dt
(3)
d.~_-- aki (l-C) dt
(4)
d_CC = ki (l-C) dt
(5)
= 3k3(1-~1) 2
(6)
dt i n -- ~1[~.~0
(7)
kl = kp(2f kiIo/kt)~
(8)
k3 = M0 kp2/kt ki = Ai exp (-El/RgT)
Initiation
(9) (10)
651 kp = Ap exp (-Ep/RgT)
Propagation
(11)
kt = At exp (-Et/RgT)
Termination
(12)
,~1, ~2 are the zeroth, first and second moments of the dead polymer. However, ~l also represents the monomer conversion, and limits can be imposed on it in optimisation that will reflect the fixed level of conversion we want to achieve. Admittedly, this model does not give a full representation of a polymerisation reaction as some assumptions have been made to reduce the complexity that comes with trying to adequately describe a polymerisation equation. Some of the assumptions are: 1. Transfer to monomer has been neglected 2. The 'gel' or Trommsdorf effect has not been taken into account. 3. Termination is by combination alone. The model equations (1) - (12) can be represented in a compact form: f(t, x'(t), x(t), u(t), v_) = O,
(13)
[to, tf]
where t is the independent variable, x_(t) ~ R n is the set of all state variables , x_'(t) denotes the derivatives of x(t) with respect to time, u_(t) E R m is a vector of control variables such as reactor temperature and coolant flow rate, and v is a vector of time invariant parameters (design variables) such as volume of the reactor. Suitable initial conditions x ( t ) are defined at time t = to. The time interval of interest is [to, tf] and the function f" R x R ~ x R ~ x R m X R p " ) R ~ is assumed to be continuously differentiable with respect to all its arguments.
3. Optimisation Problem Formulation and Constraints The optimisation problem can be described as: Given: Fixed monomer conversion and number average chain length. The temperature profile, initial concentration of initiator Optimise" The batch time for polymerisation. So as to Minimise: Subject to: Constraints on Initiator Conversion, temperature range, polydispersity values and model equations. Mathematically, the problem can be represented as, Min. tf T, Io f(t, x'(t), x_(t), u(t), v_) = 0, [to, tf] model equations s.t. m = m* Xn = X n ,
0.95
652 (100% conversion) at the end of batch time. The optimiser optimises the initial concentration needed. It is important that if possible, no trace of initiator be left in the reactor after processing so as not to discolour the finished polymer and lead to processing problems later on. By fixing these constraints, we can optimise and get a more accurate representation of a polymer product than was possible with the use of Maximum Principle (Raja, 1995).
4. Case Studies Two case studies are presented in this paper. We shall refer to them as cases A and B. In case A, we run a simulation with temperature alone as the control variable. Case B is a rerun of case A with the addition of energy balance equations.
4.1 Case study A Two sets of results are shown in Table 1 with two different values of m* and Xn*. Three control intervals are used for both cases. For comparison, in each set the results obtained by Raja are shown in italics. As can be seen, there is a significant improvement in batch time. We have also achieved closer values of m and Xn, as well as higher Initiator conversion values. The optimum Initial Initiator Concentration is also shown in Table 1. In the first case, the I0 value is actually higher but it must be said that because of the higher conversion, the amount left in the reactor at the end is less than with the maximum principle. A word must be said here about the Polydispersity (PD) values. Ideally, PD values in free radical polymerisation will be in the range 1.5 - 2.0, so our slightly greater values are well within optimal limits. Table 1 also gives the optimised temperature profiles. Figure 1 shows the monomer conversion profiles obtained by us and that obtained by Raja (1995).
4.2 Case study B In this case study, we have developed a detailed model by adding more equations to describe the interactions taking place between the reaction, the heating/cooling systems and other considerations that arise as a result of energy balance considerations. This model can be used for on-line control studies. Two of the extra equations are presented here:
d_T.T = - AH~ R__~V__ at Cp V p
= F_.j_(Tj0 -Tj) +
dt
Vj
-
O)_ CpV p
( 14 )
O_l
(15)
CpjVjpj
Coolant flow is here used as the control variable to regulate the reactor temperature. The coolant flow is bounded between 0 m3/s (no coolant) and 15 m3/s (full flow).
653 Table 1- Summary of the results of case study A. Run
m*
Xn*
i
03
500
2
IoE-2
Xn
m
499 0.9038 0.30 0.95 520 0.892 0.29 0.92 1510 1500 0.6608 0.7 0.997 1533 0.694 0.7 0.96 Optimum Temperature Profile
0.7
362.7
Temp. (K) I
365.6
Run 1 ~99.9
Time (secs). 0 Temp. (K)
Run
331.6
Po
Time(sec)
1.57
3300 4110 272336 371291
1.53 1.59 1.63 373.5
199
3300
336.7
I
343.9
I
Time (secs). 0
I
94336
I
183336
As in Run 1 of the case study A, we seek to achieve 0.3 monomer conversion and Xn of 500. The results are presented in Table 2. The optimal coolant flow profile is shown in Fig. 2. Note that in this case the batch time is about 11% higher compared to the Run 1 of Table 1. This is due to the use of a more rigorous model in the optimisation. Still, even with the addition of energy balance equations, the batch time is still lower than that obtained by Raja (1995). Table 2: Summary of the results of case study B. IoE-2 0.892
m
c 0.3
0.915
Xn
Xw
515
Fig. 1" Monomer conversion profile for Maximum Principle and technique.
805
Po 1.56
Time (secs) 3700
Fig. 2: Coolant flow profile for case study B
5. Conclusions From the results, the advantages of the use of CVP in the dynamic optimisation technique are immediately apparent. Apart from decreases in batch times and initial initiator concentration, we get closer values of m and Xn. A much lower batch time is
654 achieved using the proposed technique compared to that obtained previously by Raja (1995) using the Maximum Principle based optimisation method. This would represent cumulative savings of thousands of Euros per year for a process plant. Other advantages include: (a) the elimination of the two-step process of Raja (1995), (b) having initiator conversion close to 100% at the end of the batch which reduces the risk of fouling caused by any un-reacted initiator at the end of the reaction (c) the optimisation of the switching times (this was not possible with the Maximum Principle). In addition, from a control viewpoint, it is much easier to implement discrete timecontrol intervals than the constantly changing temperature profile obtained by Raja (1995) using the Maximum Principle. Nomenclature
a - Density of reacting mixture Ai . Pre-exponential factor for Initiator decomposition, 1.58 x 1015sec~ Ap Pre-exponential factor for Propagation, 1.1051 x 1071/g-mol-sec At . Pre-ex~ponential factor for termination, 1.255 x 10"l/g-mol-sec. C - Initiator Conversion Cp. Specific heat of reaction, 0.54 cal/g K Cpj . Specific heat of cooling liquid, 1.003 cal/g K Ei Activation Energy for Initiator decomposition, 30800 cal/g-mol Ep. Activation Energy for propagation, 7060 cal/g-mol Kp_ Rate constant for Propagation, (1/gmol- sec ) Et-Activation Energy for termination, 1680 cal/g-mol f- Initiator Efficiency, 0.6 Fj_ C o o l a n t flow rate, m3/sec
Kt. Rate constant for Termination, (l/gmol- sec ) Kn- (n=l,2..),Rate expressions for model equations. M0 Initial Monomer Concentration, 8.7006 mol/1 Rg. Universal Gas Constant, 1.987 cal/gmol.K Rp - Rate of polymerisation T. Temperature of reaction, (K) Tj_Temp. of cooling liquid, (K) Tj0 .Inlet Temp. of cooling liquid, (K) V. Volume of reactor, 300 litres Vj. Volume of Coolant, 100 litres Xn-Number Ave. Molecular Weight ~n - (n=0, 1, 2)-nth dimensionless moment for dead polymer. AHr - Enthalpy of Reaction, 16700 cal/gmol 9 - Density of reaction, 9.06 x 105g/l pj. Density of Coolant, 9.97 x 105 g/1 * - Optimal value of Indicated Parameter
I0-Initial Initiator Conc., mol/1 Ki. Rate constant for Initiation, (sec ~) References
Aziz, N and Mujtaba, I.M., 2002. Chem. Engrg. Journal, 85, 317. Fikar, M., Latifi, M.A., Fournier, F., Creff, Y., 1998. Computers Chem. Engng; 22 (supl.), pp. $ 6 2 5 - $628 Raja, T. M. 1995, Modelling and control of the Polymerisation of styrene, PhD. Thesis, University of Bradford. Ray, W.H, 1986., IEEE Control System Magazine, 6, 4, 3-8. Wallis, J.P.A., Ritter, R.A., Andre, H., 1975, Part I &II, AIChE Journal, 21, 4, 691-698. Young, R. J. and Lovell, P.A., 1991. Introduction to Polymers, 2nd edition, Chapman and Hall, N.Y.
649
Dynamic Optimisation of Styrene Polymerisation in Batch Reactors E.E. Ekpo and I.M. Mujtaba* School of Engineering Design & Technology University of Bradford, West Yorkshire BD7 1DP, UK.
Abstract This paper focuses on the optimisation of free radical polymerisation of styrene using 2,2' azobisisobutyronitrile catalyst, (AIBN) as initiator in a batch reactor. A dynamic optimisation method is used to find the optimal temperature profile that will yield a desired level of monomer conversion and number average molecular weight in minimum time. The batch time is divided into a finite number of intervals and piecewise constant temperature is used in each interval. In each interval, the temperature and the length of the interval are optimised.
Keywords:
Batch Polymerisation, Molecular Weight Distribution, Optimal Reactor Temperature, Dynamic Optimisation
1. Introduction In today's marketplace where demands can oscillate widely in short periods, versatile systems that can adapt to changing trends are needed. This is where a batch reactor has advantage over a continuous reactor. Small batch reactors require less auxiliary equipment, and their control systems are less elaborate and costly (Ray, 1986). Polystyrene is a widely used polymer with many day-to-day uses from flimsy foam packaging to more durable plastic parts used in automobiles. It is certainly in the interest of this billion-Euro-a-year business that production be optimised. In this work, a minimum time dynamic optimisation problem is formulated for free radical bulk polymerisation in a batch reactor using temperature or coolant flow as the control variable. Here, the problem is posed as a Nonlinear Programming (NLP) problem, which is solved using an SQP-based optimisation technique. The technique has been widely used in the past (Aziz and Mujtaba, 2002). The goal is to find the optimal temperature or coolant flow profiles and initial amount of initiator (I0) to yield specified values of monomer conversion (m*) and the Number Average Chain length (Xn) in minimum time. Many chemical engineering systems are transient and are usually modelled by combinations of Differential-Algebraic Equations (DAEs) of varying complexities. Suitable numerical techniques to simulate and optimise the operation of these systems have been sought in the past. Some of these are: a. The use of Dynamic Programming techniques. b. The use of NLP based Dynamic Optimisation techniques. c. The use of Pontryagin's Maximum Principle. * all correspondences to Dr. I.M. Mujtaba, email: [email protected]
650 The last two techniques are of interest to us. Raja (1995) used the method of Lagrange multipliers to estimate the initial Initiator Concentration, I0 needed and the time it would take for the polymerisation to go on to completion and then, Pontryagin's Maximum Principle to optimise the values of X, and m (a two-step process). Here, we resort to using a Control Vector Parameterisation (CVP) technique to pose the dynamic optimisation problem as NLP problem which is then solved with SQP-based optimisation routines. The main attraction of the CVP technique is its ability to handle large systems without needing to solve an excessively large optimisation problem. CVP methods are also faster than IDP methods (Fikar et.al.,1998). As a result, the need for a two-part process as in the case of the maximum principle is removed. We also get closer values for the specified parameters in shorter batch times.
2. Model Equations The model used in this work has been developed from two sources. (Wallis et al, 1975; Young and Lovell, 1991). It is a set of DAEs represented in terms of the first three moments of the dead polymer that describe the changes taking place in the batch reactor. Monomer conversion (m) and number average molecular weight (Xn) have to be evaluated from this system of DAEs. From these, the MWD can be represented by the Shulz formula: W(j) = 10 * (-log W )b+2 jb+l Tj / r(b+2)
(1)
N(j) = (W(j)/j)/Z(W(j)/(j)
(2)
where b = -(2-PD)/(1-PD) and T = exp (-(b+ 1)/Xn). W(j) and N(j) are defined as the continuous weight and number average molecular weights respectively. Some of the model equations are presented below: d__.~! = kl (l-C) ~ (1-.~!) dt
(3)
d.~_-- aki (l-C) dt
(4)
d_CC = ki (l-C) dt
(5)
= 3k3(1-~1) 2
(6)
dt i n -- ~1[~.~0
(7)
kl = kp(2f kiIo/kt)~
(8)
k3 = M0 kp2/kt ki = Ai exp (-El/RgT)
Initiation
(9) (10)
651 kp = Ap exp (-Ep/RgT)
Propagation
(11)
kt = At exp (-Et/RgT)
Termination
(12)
,~1, ~2 are the zeroth, first and second moments of the dead polymer. However, ~l also represents the monomer conversion, and limits can be imposed on it in optimisation that will reflect the fixed level of conversion we want to achieve. Admittedly, this model does not give a full representation of a polymerisation reaction as some assumptions have been made to reduce the complexity that comes with trying to adequately describe a polymerisation equation. Some of the assumptions are: 1. Transfer to monomer has been neglected 2. The 'gel' or Trommsdorf effect has not been taken into account. 3. Termination is by combination alone. The model equations (1) - (12) can be represented in a compact form: f(t, x'(t), x(t), u(t), v_) = O,
(13)
[to, tf]
where t is the independent variable, x_(t) ~ R n is the set of all state variables , x_'(t) denotes the derivatives of x(t) with respect to time, u_(t) E R m is a vector of control variables such as reactor temperature and coolant flow rate, and v is a vector of time invariant parameters (design variables) such as volume of the reactor. Suitable initial conditions x ( t ) are defined at time t = to. The time interval of interest is [to, tf] and the function f" R x R ~ x R ~ x R m X R p " ) R ~ is assumed to be continuously differentiable with respect to all its arguments.
3. Optimisation Problem Formulation and Constraints The optimisation problem can be described as: Given: Fixed monomer conversion and number average chain length. The temperature profile, initial concentration of initiator Optimise" The batch time for polymerisation. So as to Minimise: Subject to: Constraints on Initiator Conversion, temperature range, polydispersity values and model equations. Mathematically, the problem can be represented as, Min. tf T, Io f(t, x'(t), x_(t), u(t), v_) = 0, [to, tf] model equations s.t. m = m* Xn = X n ,
0.95
652 (100% conversion) at the end of batch time. The optimiser optimises the initial concentration needed. It is important that if possible, no trace of initiator be left in the reactor after processing so as not to discolour the finished polymer and lead to processing problems later on. By fixing these constraints, we can optimise and get a more accurate representation of a polymer product than was possible with the use of Maximum Principle (Raja, 1995).
4. Case Studies Two case studies are presented in this paper. We shall refer to them as cases A and B. In case A, we run a simulation with temperature alone as the control variable. Case B is a rerun of case A with the addition of energy balance equations.
4.1 Case study A Two sets of results are shown in Table 1 with two different values of m* and Xn*. Three control intervals are used for both cases. For comparison, in each set the results obtained by Raja are shown in italics. As can be seen, there is a significant improvement in batch time. We have also achieved closer values of m and Xn, as well as higher Initiator conversion values. The optimum Initial Initiator Concentration is also shown in Table 1. In the first case, the I0 value is actually higher but it must be said that because of the higher conversion, the amount left in the reactor at the end is less than with the maximum principle. A word must be said here about the Polydispersity (PD) values. Ideally, PD values in free radical polymerisation will be in the range 1.5 - 2.0, so our slightly greater values are well within optimal limits. Table 1 also gives the optimised temperature profiles. Figure 1 shows the monomer conversion profiles obtained by us and that obtained by Raja (1995).
4.2 Case study B In this case study, we have developed a detailed model by adding more equations to describe the interactions taking place between the reaction, the heating/cooling systems and other considerations that arise as a result of energy balance considerations. This model can be used for on-line control studies. Two of the extra equations are presented here:
d_T.T = - AH~ R__~V__ at Cp V p
= F_.j_(Tj0 -Tj) +
dt
Vj
-
O)_ CpV p
( 14 )
O_l
(15)
CpjVjpj
Coolant flow is here used as the control variable to regulate the reactor temperature. The coolant flow is bounded between 0 m3/s (no coolant) and 15 m3/s (full flow).
653 Table 1- Summary of the results of case study A. Run
m*
Xn*
i
03
500
2
IoE-2
Xn
m
499 0.9038 0.30 0.95 520 0.892 0.29 0.92 1510 1500 0.6608 0.7 0.997 1533 0.694 0.7 0.96 Optimum Temperature Profile
0.7
362.7
Temp. (K) I
365.6
Run 1 ~99.9
Time (secs). 0 Temp. (K)
Run
331.6
Po
Time(sec)
1.57
3300 4110 272336 371291
1.53 1.59 1.63 373.5
199
3300
336.7
I
343.9
I
Time (secs). 0
I
94336
I
183336
As in Run 1 of the case study A, we seek to achieve 0.3 monomer conversion and Xn of 500. The results are presented in Table 2. The optimal coolant flow profile is shown in Fig. 2. Note that in this case the batch time is about 11% higher compared to the Run 1 of Table 1. This is due to the use of a more rigorous model in the optimisation. Still, even with the addition of energy balance equations, the batch time is still lower than that obtained by Raja (1995). Table 2: Summary of the results of case study B. IoE-2 0.892
m
c 0.3
0.915
Xn
Xw
515
Fig. 1" Monomer conversion profile for Maximum Principle and technique.
805
Po 1.56
Time (secs) 3700
Fig. 2: Coolant flow profile for case study B
5. Conclusions From the results, the advantages of the use of CVP in the dynamic optimisation technique are immediately apparent. Apart from decreases in batch times and initial initiator concentration, we get closer values of m and Xn. A much lower batch time is
654 achieved using the proposed technique compared to that obtained previously by Raja (1995) using the Maximum Principle based optimisation method. This would represent cumulative savings of thousands of Euros per year for a process plant. Other advantages include: (a) the elimination of the two-step process of Raja (1995), (b) having initiator conversion close to 100% at the end of the batch which reduces the risk of fouling caused by any un-reacted initiator at the end of the reaction (c) the optimisation of the switching times (this was not possible with the Maximum Principle). In addition, from a control viewpoint, it is much easier to implement discrete timecontrol intervals than the constantly changing temperature profile obtained by Raja (1995) using the Maximum Principle. Nomenclature
a - Density of reacting mixture Ai . Pre-exponential factor for Initiator decomposition, 1.58 x 1015sec~ Ap Pre-exponential factor for Propagation, 1.1051 x 1071/g-mol-sec At . Pre-ex~ponential factor for termination, 1.255 x 10"l/g-mol-sec. C - Initiator Conversion Cp. Specific heat of reaction, 0.54 cal/g K Cpj . Specific heat of cooling liquid, 1.003 cal/g K Ei Activation Energy for Initiator decomposition, 30800 cal/g-mol Ep. Activation Energy for propagation, 7060 cal/g-mol Kp_ Rate constant for Propagation, (1/gmol- sec ) Et-Activation Energy for termination, 1680 cal/g-mol f- Initiator Efficiency, 0.6 Fj_ C o o l a n t flow rate, m3/sec
Kt. Rate constant for Termination, (l/gmol- sec ) Kn- (n=l,2..),Rate expressions for model equations. M0 Initial Monomer Concentration, 8.7006 mol/1 Rg. Universal Gas Constant, 1.987 cal/gmol.K Rp - Rate of polymerisation T. Temperature of reaction, (K) Tj_Temp. of cooling liquid, (K) Tj0 .Inlet Temp. of cooling liquid, (K) V. Volume of reactor, 300 litres Vj. Volume of Coolant, 100 litres Xn-Number Ave. Molecular Weight ~n - (n=0, 1, 2)-nth dimensionless moment for dead polymer. AHr - Enthalpy of Reaction, 16700 cal/gmol 9 - Density of reaction, 9.06 x 105g/l pj. Density of Coolant, 9.97 x 105 g/1 * - Optimal value of Indicated Parameter
I0-Initial Initiator Conc., mol/1 Ki. Rate constant for Initiation, (sec ~) References
Aziz, N and Mujtaba, I.M., 2002. Chem. Engrg. Journal, 85, 317. Fikar, M., Latifi, M.A., Fournier, F., Creff, Y., 1998. Computers Chem. Engng; 22 (supl.), pp. $ 6 2 5 - $628 Raja, T. M. 1995, Modelling and control of the Polymerisation of styrene, PhD. Thesis, University of Bradford. Ray, W.H, 1986., IEEE Control System Magazine, 6, 4, 3-8. Wallis, J.P.A., Ritter, R.A., Andre, H., 1975, Part I &II, AIChE Journal, 21, 4, 691-698. Young, R. J. and Lovell, P.A., 1991. Introduction to Polymers, 2nd edition, Chapman and Hall, N.Y.