Batch-Wise Nonlinear Model Predictive Control of a Gas Antisolvent

Article pubs.acs.org/IECR

Batch-Wise Nonlinear Model Predictive Control of a Gas Antisolvent Recrystallization Process for the Uniform Production of Micronized HMX with Carbon Dioxide as the Antisolvent Shin Je Lee, Sungho Kim, Bumjoon Seo, Youn-Woo Lee, and Jong Min Lee* School of Chemical and Biological Engineering, Institute of Chemical Processes, Seoul National University, Seoul 08826, Republic of Korea ABSTRACT: Novel crystallization processes that use supercritical fluids have recently attracted considerable attention because they can overcome the problems associated with conventional crystallization processes. The gas antisolvent (GAS) process is one of the promising techniques and has been applied to many applications. However, control of the GAS process is a quite challenging problem and has not yet been studied due to the complex liquid−vapor equilibrium and particle formation kinetics. This work proposes a batch-wise nonlinear model predictive control (BNMPC) approach to the GAS process to obtain the desired particle size distribution (PSD) of HMX, a widely used explosive, which should be small and uniform for stability. Although a dynamic model of the GAS crystallization is required for BNMPC, the previously developed model is too complex for real-time applications. We propose a model simplification strategy for the conventional model using experimental data. We also employ a high-resolution method (HRM) to solve effectively a partial differential equation (PDE). The simulation results show that BNMPC can produce more uniform and smaller HMX particles.

1. INTRODUCTION Crystallization plays a considerable role in producing various chemical products, such as polymers, dyes, pharmaceuticals, and explosives in fine-chemical industries. The industrial applications of crystallization include semiconductor manufacturing, drug delivery, and the food industry. Recently, a number of new applications have also involved crystallization processes, such as the production of nano and amorphous materials.1 Over the past decades, crystallization techniques have been well established in nearly all process industries as a method for producing, purifying, and recovering solid materials, and these techniques have witnessed major advancements in both academic and industrial areas. However, there are several practical problems associated with the process, which are mostly attributed to the use of large amounts of toxic solvents. The produced substances can be contaminated with the solvent, and waste solvent streams are inevitably generated in most crystallization processes. Using supercritical fluids as the solvent can overcome the drawbacks of the conventional processes because this approach uses environmentally benign solvents such as CO2. Moreover, the unique physical properties of supercritical fluids make the process easy to adjust and operable under mild conditions.2,3 In particular, the gas antisolvent (GAS) recrystallization process has attracted considerable interest because of its ability to produce fine particles. Two-way mass transfers of antisolvent and organic solvent facilitate uniform nucleation and almost instantaneous crystallization to yield ultrafine particles with a narrow particle size distribution (PSD).4 A conceptual representation of the GAS process is presented in Figure 1a, and lab-scale GAS equipment is shown in Figure 1b. A semibatch precipitator with constant volume, V, has one inlet to which the compressed CO2 gas is injected. A solution dissolved with the solute is initially loaded in a precipitator before the © 2015 American Chemical Society

Figure 1. (a) Schematic representation of the GAS process and (b) experimental equipment.

process starts. The antisolvent CO2 is pumped to high pressure over its critical point and injected into the vessel from the bottom to achieve better mixing of the solvent and antisolvent. After the high-pressure CO2 is injected, the solution volume is expanded; however, the precipitator volume is maintained constant for the entire process. Acetone does not vaporize, and thus, it does not affect the supersaturation or particle formation. When the CO2 dissolves in acetone, the solvent strength on the solute is reduced, the system becomes supersaturated and particles are produced. After a holding time, the solution is drained under isobaric conditions to wash and collect the particles. With increasing interest, many studies on the GAS process have been reported. Thermodynamic studies of phase behavior Received: Revised: Accepted: Published: 11894

May 6, 2015 November 2, 2015 November 5, 2015 November 6, 2015 DOI: 10.1021/acs.iecr.5b01690 Ind. Eng. Chem. Res. 2015, 54, 11894−11902

Article

Industrial & Engineering Chemistry Research

algebraic equations, and it is very complex to solve explicitly at every sample time. Therefore, we propose a model simplification approach by simple algebraic manipulations and using experimental data. Then, we use a high-resolution method (HRM) as the method for solving the PDE in the GAS model because it is more stable and accurate for a steep function than FDM or FEM. The remainder of this paper is organized as follows. In Section 2, a fundamental model of the GAS process and its simplification approach using the experimental data are presented, and the HRM discretization is also described. The BNMPC algorithm is presented in Section 3, and the simulation results of HMX recrystallization are shown in Section 4, followed by concluding remarks in Section 5.

in the GAS process and their experimental results have been explored.5−7 Some studies have examined the accurate determination of the volume expansion in the GAS process because it is related to the product quality.4,8,9 Many experimental studies on selecting the optimal operating conditions have been performed, and key parameters include the antisolvent addition rate, temperature, organic solvent type, and initial solution concentration.10−14 Systematic mathematical modeling of the GAS process was developed by Muhrer et al. 15 using a population balance model (PBM) with thermodynamics and particle formation kinetics. The developed model has been validated with experimental data in the literature,16,17 and a novel discretization method for solving the model has also been reported.17 This work addresses the recrystallization of cyclotetramethylenetetranitramine (HMX), a widely used explosive, using the GAS process with acetone and CO2 as a solvent and antisolvent, respectively. It has already been shown that HMX can be recrystallized into small and uniform particles using the GAS process in the literature.18,19 High-energy materials such as explosives must be small with a uniform PSD for performance and stability. This emphasizes the need for controlling the GAS process. Although obtaining fine crystals with a uniform distribution is a critical design problem in the GAS process, studies on controlling the GAS process have hardly progressed compared with the extensive studies for model development and various applications. This is because controlling the PSD of the GAS process is a challenging problem due to its inclusion of complex vapor−liquid equilibrium and particle kinetics. Although there are a few studies on the control of liquid antisolvent crystallization,20−22 the liquid antisolvent process is quite different from the GAS process in that it does not include vapor−liquid equilibrium in the model. This work proposes a nonlinear model predictive control (NMPC) approach for controlling the PSD of the GAS process because the system model is highly nonlinear. MPC applications for conventional crystallization have been widely studied and implemented in industries.23−25 In industrial crystallization processes, the control problem has always been an issue and has significantly advanced due to the importance of manufacturing crystals with consistent quality. These advances have been enabled by the progress in in situ realtime sensor technologies.26 Robust performance can be achieved by repeating the optimization online on the basis of real-time measurements and state estimation.27 The nonlinear GAS model is linearized every sample time in the control algorithm, and the conventional MPC is implemented batchwise because the GAS recrystallization is a batch process. Simulation results show that the final product quality of HMX is improved at the end of the final batch by a batch-wise nonlinear model predictive control (BNMPC) strategy. To the authors’ best knowledge, there are no published results regarding the implementation of MPC to the GAS recrystallization, which is more challenging because of the complexity raised by the phase equilibrium, crystallization kinetics, and calculation of a partial differential equation (PDE) describing population balances in the model. In the MPC algorithm, a dynamic process model of the system is used to predict the effects of future actions on the outputs. Therefore, the process model is the essential element of the MPC controller. However, the previously developed model of the GAS process consists of a PDE, ODEs, and

2. DYNAMIC MODEL OF THE GAS PROCESS AND MODEL SIMPLIFICATION A mathematical model for the GAS process was completely developed by Muhrer et al.15 using a population balance model (PBM) combined with particle formation dynamics and thermodynamics. After introducing this model, we propose a simplification strategy of the model to be effectively applied to control schemes. Then, a high-resolution method (HRM) is proposed to solve the PBM because the HRM exhibits superior numerical accuracy compared to other well-known discretization methods such as FDM or FEM. 2.1. Fundamental Model of the GAS Process. The PBM is used to express how particle sizes are distributed as a partial differential equation (PDE). The one independent variable is the time, t; the other variable is a property coordinate, such as the particle size, L. General crystallization dynamics are described by28 ∂{G(L , t )n(L , t )} ∂n(L , t ) = q(L , t , f ) + ∂L ∂t

(1)

where n(L,t) is the population density function representing a particle size distribution (PSD); t and L are the time and the size, respectively; G(L,t) is the growth rate of particles; and q(L,t,f) is the creation/depletion rate, which accounts for nucleation, aggregation, breakage, and material leaving or entering the system. In the dynamic GAS model based on eq 1, several assumptions are made while retaining the basic dynamic behavior of the system:15 the temperature in the vessel is maintained constant such that no energy balance is needed; the growth rate, G, is size independent; and aggregation and breakage contributions of particles are also neglected. These assumptions simplify eq 1 to ∂n ∂n +G =0 ∂t ∂L

(2)

The change of the PSD according to the volumetric expansion of the liquid phase should be added to eq 2 because it determines the supersaturation of the solute, which is the driving force for particle formation ∂n ∂n n d(NLvL) +G + =0 ∂t ∂L NLvL dt

(3)

where NL [mol] is the molar hold-up in the liquid phase and vL [m3/mol] is the molar volume of the liquid phase. Therefore, NLvL is the liquid volume, and d(NLvL)/dt represents the liquid volume change with time. The last term of eq 3 explains the 11895

DOI: 10.1021/acs.iecr.5b01690 Ind. Eng. Chem. Res. 2015, 54, 11894−11902

Article


well describes the vapor−liquid equilibrium of the solvent and antisolvent with a moderate complexity in the GAS process. Twelve unknown variables (n, NL, NV, NP, vL, vV, xA, xS, xP, yA, yS, and P) can be obtained by solving a system of partial differential, ordinary differential, and algebraic equations of eqs 3−11 and 13−15 with the following initial conditions

effect of the liquid volume expansion on the PSD in the GAS process. Finally, the PBM for the GAS process can be obtained as the PDE. The material balances on the antisolvent (A), solvent (S), and particle (P) are d(NLxA + NVyA )

= QA

dt d(NLxS + NVyS ) dt

=0

d(NLx P + NP) =0 dt

(4)

(5)

(6)

NLvLk vm3 vP

yA + yS = 1

(9)

(10)

fS,L = fS,V

(11)

(i = A, S, α = L, V)

(15)

(21)

⎧ B′ + B ′′ S > 1 B=⎨ ⎩0 S≤1

(22)

⎧ kg(S − 1)g S > 1 G=⎨ ⎩0 S≤1

(23)

fP,L fP,P

(24)

⎡ v (P − P0) ⎤ fP,P = fP,L (P0 , T , x0) × exp⎢ P ⎥ ⎣ ⎦ RT

where the liquid and vapor molar volumes, vL and vV [m3/mol], respectively, are calculated from an equation of state as

P = P(v V , T , y)

B G

The fugacity of the solid at the solid phase, f P,P, is calculated using the Poynting correction factor

(13)

(14)

(20)

S=

(12)

P = P(vL , T , x)

x PNL = NP0

where B′ and B″ are the primary and secondary nucleation rates, respectively; S is the supersaturation level; and kg and g are the rate constants. Constitutive equations of the primary and secondary nucleation rates can be found in the literature.15,30 The nucleation and growth rates are functions of S, and the control of crystallization is to balance the nucleation and growth rates to adjust the particle size and achieve the desired PSD. The supersaturation, S, is defined by the ratio of fugacities of the solid in the liquid and solid phases,

where ϕ is the fugacity coefficient. The total volume, V, of the crystallizer is constant and satisfies NLvL + NVv V = V

(19)

⎪

The fugacities in the liquid and vapor phases are calculated as fi , α = zi , αϕi , αP

xsNL + yS NV = NS0

⎪

and liquid−vapor equilibrium leads to the isofugacity relationships fA,L = fA,V

(18)

where B and G are the nucleation and growth rates, respectively, defined as

where kv, vP, and m3 are the volume shape factor, molar volume of the solid phase, and third moment of the population density function, respectively. The mole fractions of the three components are subject to (8)

n(0, L) = 0

n(t , 0) =

(7)

xA + xS + x P = 1

(17)

where N0S and N0P are the initial molar amounts of solvent and solute, respectively. The boundary condition for the PDE is given as

where NV [mol] and NP [mol] are the molar hold-ups in the vapor and solid phases, respectively; xi and yi are the mole fractions of component i in the liquid and vapor phases, respectively (i = A, S, P); and QA [mol/s] is the molar flow rate of the antisolvent and is the input variable. It is found that the antisolvent addition rate, QA, has the strongest influence on the final product,15 and it can be manipulated to control particle size, PSD, and morphology.29 The molar hold-up in the solid phase, NP, is given by NP =

P = Patm

(25)

where P0 and x0 indicate a reference pressure and composition at the reference pressure, respectively. 2.2. Control-Relevant GAS Model. The comprehensive GAS model is difficult to solve explicitly because complex particle dynamics and thermodynamics are involved in the system. We present a simplified model with some algebraic manipulations by using experimental data. The model should be simplified to be solved within a sample time in the modelbased control algorithm because solving the entire mathematical model of GAS is very time-consuming. After the model is simplified, it reduces to six equations and six unknown variables when the input is fixed. Summation of eqs 4−6 using eqs 8−9 provides

The Peng−Robinson equation of state (PR-EOS) with quadratic mixing rules is used aα RT (α = L, V) P= − 2 vα − bα vα + 2vαbα − bα2 (16) where R is the gas constant, T is the temperature, vα is the molar volume of the α phase, and aα and bα are the parameters of the PR-EOS of the α phase. The PR-EOS is used because it 11896


Article

Industrial & Engineering Chemistry Research dNL d NV dNP + + = QA dt dt dt

the PBM is particularly challenging because the population density function, n(L,t), extends over orders of magnitude, and thus, the distribution is very sharp. We use the HRM because it provides high-order accuracy with a lower computational load than other finite difference or finite element methods (FDM or FEM) for sharp distributions.31,32 It does not create numerical diffusion or dispersion as in first-order and second-order methods.31,32 The HRM with second-order accuracy discretizes eq 2 to32

(26)

and we specify three derivatives as f, g, and h as follows

dNL =f dt

(27)

d NV =g dt

(28)

dNP =h dt

dni G (n L + − n Li−) = 0, + dt ΔL i

(29)

(30)

n Li+ = ni +

Additionally, the derivatives of vL and vV are assigned to p and q, respectively (31)

dv V =q dt

(32)

θi+ =

ni − ni − 1 + ϵ ni + 1 − ni + ϵ

(40)

θi− =

ni − 1 − ni − 2 + ϵ ni − ni − 1 + ϵ

(41)

ϕ(r ) = (33)

|r | + r 1 + |r |

(42)

Thus, eq 36 is changed into the discretized form of ∂ni G n =− (n Li+ − n Li−) − (vLf + NLp) ∂t ΔL NLvL

(34)

(43)

The dimension of the equation is dependent on the mesh size, i. The HRM is compared with the first-order FDM for eq 2 with the size-independent growth rate (G = 1.0 μm/s), and the boundary condition is given as n(0,t) = 0. The initial distribution is33

where h can be formulated from eq 7. The liquid phase change part of eq 3 can be separately expressed as dv ⎞ n ⎛ dNL ∂n ∂n ⎜v L +G + + NL L ⎟ = 0 ⎝ NLvL dt dt ⎠ ∂t ∂L

(39)

where ϵ = 10−10. Here, we choose the Van Leer flux limiter because it does not exhibit any numerical dispersion for onedimensional problems and provides full second-order accuracy32

and substitution of eqs 27, 28, 31, and 32 into eq 33 and rearrangement for f provides 1 (NLp + NVq + v V (Q A − h)) f=− vL − v V

1 ϕ(θi−)(ni − ni − 1) 2

(38)

and ΔL is the mesh element width; N is the divided mesh number; and ϕ(r) is the flux limiter function. It depends on the degree of smoothness of the distribution, which is defined by

We can calculate p and q from the experimental results and the fact that the liquid and vapor mole fractions, xi and yi, depend on the antisolvent addition rate, QA.5 It is assumed that the pressure change is linear with the antisolvent addition rate, and this assumption is found to be reasonable by solving the equations describing the relationship between the pressure change and addition rate. Then, the molar volumes of the liquid and vapor phases, vL and vV, can be calculated from eq 16 with mole fractions, and the pressure and derivatives are also obtained. The differentiation of eq 13 leads to dvL dN dv dN + v L L + NV V + v V V = 0 dt dt dt dt

1 ϕ(θi+)(ni + 1 − ni) 2

n Li− = ni − 1 +

dvL =p dt

NL

(37)

where

Then, eq 26 can be represented as f + g + h = QA

i = 1, ..., N

(35)

and substitution of eqs 27 and 31 into eq 35 provides

⎧1 × 1010 if 10 μm < L < 20 μm n(L , 0) = n0(L) = ⎨ ⎩0 otherwise ⎪

∂n ∂n n = −G − (vLf + NLp) ∂t ∂L NLvL

⎪

(36)

(44)

Finally, the simplified model consists of eqs 27−29, 31, 32, and 36 with six unknown variables of NL, NV, NP, vL, vV, and n. The initial and boundary conditions are the same as those of the fundamental model. Then, the partial derivative of n with respect to L in eq 36 should be discretized to solve the PDE, and we apply a high-resolution method (HRM) to the PDE in the next subsection. 2.3. HRM Solution Strategy. A specific numerical solution scheme is required to solve numerically the GAS model because it includes the PDE in the PBM. The numerical simulation of

The analytical solution of this problem with an initial profile is the initial profile translated by a length Gt, that is, n(L , t ) = n0(L − Gt )

(45)

The population densities for three solution approaches are compared in Figure 2. The HRM exhibits better accuracy than the FDM, and it does not considerably change as the simulation time increases, whereas the accuracy of FDM decreases as time increases due to the numerical dispersion. 11897


Article


which can be used to compute the control actions as in the usual linear MPC formulations, is developed using the following equations.34 For a prediction horizon p, Ykj + 1 | k = Y k0,+j 1 | k + Skj ΔUkj , j = 1, ..., p

(54)

where

Ykj + 1 | k Figure 2. Comparison of HRM and FDM for the solution of eq 45 at (a) 30 s and (b) 60 s.

(46)

yk = g (xk , uk) + vk

(47)

and Sjk is the matrix of step response coefficients for the linear system, which may be calculated as described in the literature.35 ykj + l|k is the predicted value of y at time k + l based on the information available at time k. The vector Y0,j k+1|k provides the output trajectories if no control actions were taken. 3.3. Calculation of Control Action. At time k, the control action to be solved online at every sampling instance in the NMPC algorithm is

where xk, uk, and yk are the vectors of nx system states, nu inputs, and ny measured variables at the kth sampling instance and wk and vk are the zero-mean white noises on the system states and measured variables. The system dynamics and the measurements are described by the nonlinear functions f and g. The states, inputs, and outputs of the GAS recrystallization process are x = [ni , NL , NV , NP , vL , v V ]T , i = 1, ..., M

minΔU{||Q (R kj + 1 | k − Ykj + 1 | k)||22 + ||R ΔUkj||22 } j j umin ≤ ukj+ l | k ≤ umax , j Δukj+ l | k ≤ Δumax , j j ymin ≤ ykj+ n | k ≤ ymax ,

(49)

∑ niΔL

where M is the mesh size of the particle density function, n. Therefore, the dimension of the system is (M+5) . We set M as 100 and ΔL as 1 μm. The measured output, y, represents the total number of particles with particle sizes from 0 to 100 μm. The calculation of the total number of particles is from the moment of a distribution. At time k, the nonlinear functions f and g in xk and uk are linearized with respect to xk and uk−1 using a Jacobian matrix at every sample time. xk̅ + 1 = Ak xk̅ + Bk uk̅ (51)

yk̅ = Ckxk̅

u0j + 1 = uNj , j = 1, ..., J

(52)

∂f |x , u , ∂x k k −1

Bk =

∂f |x , u , ∂u k k −1

Ck =

n = 1, ..., p

(57)

(58)

where J is the total batch number. Therefore, the PSD quality approaches the set point trajectory as the number of batch processes increases.

where Ak =

l = 0, ..., m − 1

where is the vector of future set points corresponding to Yjk+1|k in the jth batch, Q and R are the weighting matrices, and umin, umax , ymin, and ymax are the minimum and maximum values of the input and output variables, respectively. User-selected tuning parameters are the prediction horizon (p), control horizon (m)(m ≤ p), and weighting matrices in the objective function (Q and R). The above constrained optimization problem can be expressed as a quadratic program (QP), which is solved online at every sampling time. The first element of ΔU is implemented on the plant, and the optimization problem is reformulated and solved at the next time step iteratively utilizing the most current measurements. 3.4. Batch Input Update. After the jth batch process is finished, the final input values become the initial input values in the next (j+1)th batch as follows

(50)

i=1

l = 0, ..., m − 1

Rjk+1|k

M

y=

(56)

subject to

(48)

u = QA

⎤ ⎡ Δu j = u j − u j k|k k|k k−1|k ⎥ ⎢ ⎢ Δu j = ukj+ 1 | k − ukj| k ⎥ ⎥ ΔUkj = ⎢ k + 1 | k ⎥ ⎢ ⋮ ⎥ ⎢ ⎥⎦ ⎢⎣ Δukj+ m − 1 | k (55)

3. BATCH-WISE NONLINEAR MODEL PREDICTIVE CONTROL (BNMPC) 3.1. Nonlinear Process Model and Linearization. The nonlinear process model has the following representation xk + 1 = f (xk , uk) + wk

⎡yj ⎤ ⎢ k+1|k ⎥ ⎢ j ⎥ y = ⎢ k + 2 | k ⎥, ⎢ ⋮ ⎥ ⎢ ⎥ ⎢yj ⎥ ⎣ k+p|k ⎦

∂g |x , u ∂x k k −1

4. RESULTS AND DISCUSSION 4.1. Open-Loop Simulation and Experimental Validation. The dynamic model of the GAS process is simulated for 100 s with constant CO2 addition rates of 20, 50, and 100 mL/min to determine the effect of the addition rate on the final PSD. The temperature is assumed to be fixed at 30 °C because the GAS process is an isothermal process. The particle size, L, is discretized from 0 to 100 μm with a mesh size of 1 μm, and thus, the total number of meshes is 101. Calculations are

(53)

and x¯k = xk − xs, u¯k = uk − us, and y¯k = yk − ys. After linearization of nonlinear functions, the conventional MPC algorithm can be applied to the linearized state-space model. 3.2. Linear Prediction of Future Outputs. Because the system is a semibatch process, the MPC algorithm is implemented within each batch. This is possible by using the real-time measurements of PSD using ATR-FTIR or ATR-UV/ V.26 In the jth batch, a linear prediction of future outputs, 11898


Article

Industrial & Engineering Chemistry Research performed using MATLAB R2014a with an Intel(R) Core(TM) i5-3570 CPU @ 3.40 GHz (4 CPUs). Figure 3 shows the

Figure 3. Final PSDs of HMX at 20, 50, and 100 mL/min at 30 °C.

open-loop simulation results of the final PSDs when the CO2 addition rates are 20, 50, and 100 mL/min. It is observed that the PSD becomes narrower and that the average particle size decreases as the CO2 addition rate increases. This result is because the fast CO2 addition rate induces a high supersaturation level in a short time, leading to a sudden burst of nucleation while simultaneously preventing the particle growth.19 In other words, the fast CO2 addition rate increases the nucleation-to-growth ratio, and this ratio varies with the antisolvent addition rates. Therefore, with the rapid CO2 addition rate, a large amount of small particles that do not achieve full growth are formed at once. The mean particle sizes and variances of the PSDs are also provided in Table 1 to compare quantitatively the results. From the variances, we can observe that the PSD becomes narrower and has uniform particles as the CO2 addition rate increases.

Figure 4. (a−d) FESEM of HMX at 25, 30, 40, and 50 °C and CO2 addition rate of 10 mL/min and (e−h) FESEM of HMX at 25, 30, 40, and 50 °C and CO2 addition rate of 50 mL/min.

Table 1. Mean Particle Sizes and Variances of the Final PSDs CO2 addition rate (mL/min)

mean size (μm)

variance

20 50 100

44.0 38.9 28.5

501 342 106

addition rate is 50 mL/min, the temperature has a strong influence on the particle size, as shown in Figure 5b, whereas under 10 mL/min, the particle size does not considerably change, as shown in Figure 5a. Moreover, when the operating temperature is greater than 40 °C, the mean particle size does not considerably decrease, even when the CO2 addition rate increases from 10 to 50 mL/min. From the open-loop test and the experimental results, it is observed that uniform particles with improved quality can be obtained by controlling the process. 4.2. Simulation Results of BNMPC. The proposed BNMPC is applied to the GAS crystallization process to obtain the desired PSD of HMX. The tuning parameters and process constraints for the BNMPC are presented in Table 2. The minimum input value is selected as 0.0134 mL/min such that the bimodal PSD is avoided, and the maximum is an arbitrary realistic value. The simulation time and sampling time are 100 and 0.5 s, respectively. Figure 6 shows the mean particle sizes and input profiles based on the results of BNMPC of the GAS process during 5 batches. In the first batch, the objective is to maximize the particle number as close to 109 as possible to make the particles small because the particle size would be smaller if a large number of particles are formed. This strategy is based on the fact that we can only measure the total particle number in the sample time.

The simulation results are also compared with the experimental data to show that the simplification of the dynamic model of GAS is acceptable. In the experiments, the CO2 addition rate is set to 10 or 50 mL/min at constant temperatures of 25, 30, 40, and 50 °C and γ-butyrolactone (GBL) is used as a solvent. Figure 4 shows the FESEM (fieldemission scanning electron microscopy) results of HMX. From Figure 4, we can observe that the real shape of HMX is not a circle, although we assume that the shape of HMX is circular in the model, and small and uniform particles can be obtained under 25 °C and 50 mL/min, as shown in Figure 4e. In Figure 5, the comparison results of mean particle sizes between the experiment and the model are shown. The observation that the operating condition of 25 °C and 50 mL/min produces a smaller and more uniform particle can be quantitatively confirmed in Figure 5. Under this condition, the mean particle size is 20 μm. We can also determine the effect of temperature on the particle size from the experiment. When the CO2 11899


Article


batch and so forth. The input profiles during each batch are slightly changed, and the final mean particle size is reduced through 5 batches. As shown in Figure 6, the input profiles suddenly increase for the first 3 batches. This result is because the particle becomes larger with time and the CO2 addition rate should increase to control the particle size to be as small as possible. Under the initial conditions, the particle size is very small, and as the particle grows, the mean particle size gradually increases. After the particle grows to some extent, the input rate suddenly increases in the middle of the process to control the particle size as small. After 3 batches, the CO2 addition rate does not considerably increase because several initial batches have the greatest effect on the PSD. The final particle sizes are quantitatively compared in Figure 7, along with the initial input values for 5 batches. The final

Figure 5. Mean particle size comparison of experimental data and simulation results at CO2 addition rates of (a) 10 and (b) 50 mL/min.

Table 2. Tuning Parameters and Constraints for the BNMPC parameter

value

p m Q R umin umax Δumax

5 5 0.3 0.2 0.0134 mL/min 1340 mL/min 0.267 mL/min

Figure 7. Comparison of final mean particle sizes and initial input values for 5 batches.

mean particle sizes for 5 batches are 33.6, 31.0, 24.7, 18.8, and 19.5 μm, and the mean size no longer decreases after the fourth batch. Additionally, the initial input values are 10.01, 0.0697, 0.0290, 0.0211, and 0.0194 mL/min. The final PSDs at 5 batches are shown in Figure 8 to show the uniformity of the PSD of HMX in the fifth batch. From Figure 8, it is observed that the final PSD at the fifth batch is the most uniform, while its mean size, 19.5 μm, is somewhat larger than 18.8 μm, the mean size of the fourth batch. In

The initial CO2 addition rate of the first batch is 10 mL/min, and it decreases to 0.0697 mL/min at the end of the batch. The final input value is maintained as the initial value of the second

Figure 6. Mean particle sizes and input profiles obtained from BNMPC for 5 batches. 11900


Article


considerable influence on the final PSD. The average particle sizes of Seeds 1, 2, and 3 are 5, 10, and 20 μm, respectively, and they are changed to 38.4, 46.1, and 53.9 μm at the end of the batch through the GAS process.

5. CONCLUDING REMARKS Control issues of the GAS process are quite challenging because the system is highly nonlinear and includes complex liquid−vapor equilibrium and particle kinetics, and it has not yet been studied. In this work, we propose a batch-wise nonlinear model predictive control (BNMPC) to control the GAS recrystallization process to make HMX, which is widely used as an explosive, small and uniform for stability. To implement the MPC, a dynamic model of the GAS process is required. The dynamic GAS model can be represented as a PBM and includes a PDE that requires a discretization method to be numerically solved. HRM is introduced because it provides high-order accuracy. Additionally, in the GAS model, several constitutive equations are also included: thermodynamics related to liquid volume expansion and nucleation and growth rates of particle kinetics. The full dynamic model of the GAS process is difficult to solve in the sampling time. Therefore, we propose a model simplification strategy using experimental data and calculation of the PR-EOS. The particle size calculated from the simplified model shows insignificant error with the experimental results. Through an open-loop simulation and the experimental results, we can observe that the particle size of HMX decreases and that its PSD is more uniform as the CO2 addition rate increases. The BNMPC algorithm is presented to control the GAS process. Within a batch, the conventional MPC is implemented, and the final input value is transferred to the initial input of the next batch. The BNMPC results show that the final mean particle size of the HMX decreases after 5 batches and that the final PSD at the fifth batch is the most uniform. In this work, we apply the MPC algorithm to the GAS crystallization process as the first step of model-based control of this system. A batch-to-batch control strategy can also be introduced to the GAS process in future research.

Figure 8. Final PSDs of HMX for 5 batches.

conclusion, fine and uniform HMX particles can be obtained from the BNMPC algorithm in the GAS crystallization process. We also observe that the PSD of the fifth batch with an input value of approximately 25 mL/min is considerably smaller than the PSD with an input value of 100 mL/min in the open-loop simulation. This result is because the open-loop process is just a one-batch process with an input value of 100 mL/min, whereas the closed-loop process goes through five batch processes, each with a specific input profile. Repeating 5 batches should be more effective than simply one batch even though it is with the maximum input rate, 100 mL/min. If we compare the PSD with 100 mL/min in Figure 3 and the PSD of the first batch in Figure 8, then the 100 mL/min clearly provides substantially smaller particles than the input profile of the first batch in the closed-loop simulation. However, through five batch processes, the PSDs are effectively controlled to provide an optimal PSD, as shown in Figure 8. After several batches, the particle sizes become smaller, and they are used as the initial particles in the next batch. Although they are dissolved in the solvent, it affects the initial PSD of HMX, as presented by Kim et al.18 In other words, the BNMPC control results are affected by the CO2 feed rate and by the initial seed PSD in a batch. To determine the effect of the initial seed PSD on the final PSD, three HMX seed PSDs with different shapes are used while the CO2 addition rate remains at a constant value of 50 mL/min. The PSD variation at the end of the batch is shown in Figure 9. It is observed that the shape of the initial PSD has a

Figure 9. Final PSDs of HMX for 5 batches. 11901


Article


■

dipropionate using carbon dioxide. Ind. Eng. Chem. Res. 2007, 46, 8009−8017. (17) Bakhbakhi, Y. A discretized population balance for particle formation from gas antisolvent process: The combined Lax-Wendroff and Crank-Nicholson method. Comput. Chem. Eng. 2009, 33, 1132− 1140. (18) Kim, S.-J.; Lee, B.-M.; Lee, B.-C.; Kim, H.-S.; Kim, H.; Lee, Y.W. Recrystallization of cyclotetramethylenetetranitramine (HMX) using gas anti-solvent (GAS) process. J. Supercrit. Fluids 2011, 59, 108−116. (19) Lee, B.-M.; Kim, S.-J.; Lee, B.-C.; Kim, H.-S.; Kim, H.; Lee, Y.W. Preparation of micronized β-HMX using supercritical carbon dioxide as antisolvent. Ind. Eng. Chem. Res. 2011, 50, 9107−9115. (20) Sheikhzadeh, M.; Trifkovic, M.; Rohani, S. Real-time optimal control of an anti-solvent isothermal semi-batch crystallization process. Chem. Eng. Sci. 2008, 63, 829−839. (21) Nagy, Z. K.; Fujiwara, M.; Braatz, R. D. Modelling and control of combined cooling and antisolvent crystallization processes. J. Process Control 2008, 18, 856−864. (22) Nowee, S. M.; Abbas, A.; Romagnoli, J. A. Model-based optimal strategies for controlling particle size in antisolvent crystallization operations. Cryst. Growth Des. 2008, 8, 2698−2706. (23) Rohani, S.; Haeri, M.; Wood, H. Modeling and control of a continuous crystallization process Part 2. Model predictive control. Comput. Chem. Eng. 1999, 23, 279−286. (24) Dokucu, M. T.; Park, M.-J.; Doyle, F. J., III Multi-rate model predictive control of particle size distribution in a semibatch emulsion copolymerization reactor. J. Process Control 2008, 18, 105−120. (25) Hermanto, M. W.; Chiu, M.-S.; Braatz, R. D. Nonlinear model predictive control for the polymorphic transformation of L-glutamic acid crystals. AIChE J. 2009, 55, 2631−2645. (26) Nagy, Z.; Braatz, R. Advances and New Directions in Crystallization Control. Annu. Rev. Chem. Biomol. Eng. 2012, 3, 55−75. (27) Larsen, P.; Rawlings, J.; Ferrier, N. An algorithm for analyzing noisy, in situ images of high-aspect-ratio crystals to monitor particle size distribution. Chem. Eng. Sci. 2006, 61, 5236−5248. (28) Ramkrishna, D. Population balances: Theory and applications to particulate systems in engineering; Academic Press: London, 2000. (29) Gallagher, P. M.; Coffey, M.; Krukonis, V.; Hillstrom, W. Gas anti-solvent recrystallization of RDX: formation of ultra-fine particles of a difficult-to-comminute explosive. J. Supercrit. Fluids 1992, 5, 130− 142. (30) Mersmann, A. Crystallization technology handbook. Drying Technol. 1995, 13, 1037−1038. (31) Gunawan, R.; Fusman, I.; Braatz, R. High resolution algorithms for multidimensional population balance equations. AIChE J. 2004, 50, 2738−2749. (32) Qamar, S.; Elsner, M.; Angelov, I.; Warnecke, G.; SeidelMorgenstern, A. A comparative study of high resolution schemes for solving population balances in crystallization. Comput. Chem. Eng. 2006, 30, 1119−1131. (33) Qamar, S.; Ashfaq, A.; Warnecke, G.; Angelov, I.; Elsner, M.; Seidel-Morgenstern, A. Adaptive high-resolution schemes for multidimensional population balances in crystallization processes. Comput. Chem. Eng. 2007, 31, 1296−1311. (34) Ricker, N.; Lee, J. Nonlinear model predictive control of the Tennessee Eastman challenge process. Comput. Chem. Eng. 1995, 19, 961−981. (35) Garcia, C.; Prett, D.; Morari, M. Model predictive control: theory and practice-a survey. Automatica 1989, 25, 335−348.

AUTHOR INFORMATION

Corresponding Author

*J. M. Lee. Tel.: +82 2 8801878. Fax: +82 2 8887295. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The Energy Efficiency & Resources Core Technology Program of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) granted financial resources from the Ministry of Trade, Industry & Energy, Republic of Korea (2012T100201687). This work is supported by the Korea Ministry of Environment as a Project for Developing EcoInnovation Technologies (GT-11-G-02-001-3).

■

REFERENCES

(1) Larsen, P.; Patience, D.; Rawlings, J. Industrial crystallization process control. IEEE Control Syst. Mag. 2006, 26, 70−80. (2) Jung, J.; Perrut, M. Particle design using supercritical fluids: literature and patent survey. J. Supercrit. Fluids 2001, 20, 179−219. (3) Knez, Z.; Weidner, E. Particles formation and particle design using supercritical fluids. Curr. Opin. Solid State Mater. Sci. 2003, 7, 353−361. (4) Mukhopadhyay, M. Partial molar volume reduction of solvent for solute crystallization using carbon dioxide as antisolvent. J. Supercrit. Fluids 2003, 25, 213−223. (5) Dixon, D. J.; Johnston, K. P. Molecular thermodynamics of solubilities in gas antisolvent crystallization. AIChE J. 1991, 37, 1441− 1449. (6) Shariati, A.; Peters, C. Measurements and modeling of the phase behavior of ternary systems of interest for the GAS process: I. The system carbon dioxide+ 1-propanol+ salicylic acid. J. Supercrit. Fluids 2002, 23, 195−208. (7) de la Fuente, J. C.; Shariati, A.; Peters, C. J. On the selection of optimum thermodynamic conditions for the GAS process. J. Supercrit. Fluids 2004, 32, 55−61. (8) De la Fuente Badilla, J.; Peters, C.; de Swaan Arons, J. Volume expansion in relation to the gas-antisolvent process. J. Supercrit. Fluids 2000, 17, 13−23. (9) Elvassore, N.; Bertucco, A.; Di Noto, V. On-line monitoring of volume expansion in gas-antisolvent processes by UV-vis spectroscopy. J. Chem. Eng. Data 2002, 47, 223−227. (10) Müller, M.; Meier, U.; Kessler, A.; Mazzotti, M. Experimental study of the effect of process parameters in the recrystallization of an organic compound using compressed carbon dioxide as antisolvent. Ind. Eng. Chem. Res. 2000, 39, 2260−2268. (11) Muhrer, G.; Mazzotti, M.; Müller, M. Gas antisolvent recrystallization of an organic compound Tailoring product PSD and scaling-up. J. Supercrit. Fluids 2003, 27, 195−203. (12) Kim, C.-K.; Lee, B.-C.; Lee, Y.-W.; Kim, H. S. Solvent effect on particle morphology in recrystallization of HMX (cyclotetramethylenetetranitramine) using supercritical carbon dioxide as antisolvent. Korean J. Chem. Eng. 2009, 26, 1125−1129. (13) Lee, B.-M.; Jeong, J.-S.; Lee, Y.-H.; Lee, B.-C.; Kim, H.-S.; Kim, H.; Lee, Y.-W. Supercritical antisolvent micronization of cyclotrimethylenetrinitramin: Influence of the organic solvent. Ind. Eng. Chem. Res. 2009, 48, 11162−11167. (14) Widjojokusumo, E.; Veriansyah, B.; Youn, Y.-S.; Lee, Y.-W.; Tjandrawinata, R. R. Co-precipitation of loperamide hydrochloride and polyethylene glycol using aerosol solvent extraction system. Korean J. Chem. Eng. 2013, 30, 1797−1803. (15) Muhrer, G.; Lin, C.; Mazzotti, M. Modeling the gas antisolvent recrystallization process. Ind. Eng. Chem. Res. 2002, 41, 3566−3579. (16) Dodds, S.; Wood, J. A.; Charpentier, P. A. Modeling of the gasantisolvent (GAS) process for crystallization of beclomethasone 11902


Batch-Wise Nonlinear Model Predictive Control of a Gas Antisolvent

Recommend Documents