Dynamic Modeling and Optimization of Batch Crystallization of Sugar

Article pubs.acs.org/IECR

Dynamic Modeling and Optimization of Batch Crystallization of Sugar Cane under Uncertainty Eusebio Bolaños-Reynoso,*,† Kelvyn B. Sánchez-Sánchez,† Galo R. Urrea-García,† and Luis Ricardez-Sandoval‡ †

División de Estudios de Posgrado e Investigación, Instituto Tecnológico de Orizaba, Avenida Oriente 9 Núm. 852 Colonia Emiliano Zapata, Orizaba, Veracruz 94320, México ‡ Department of Chemical Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada ABSTRACT: This work presents a study on the agitation rate effects on the average diameter (% volume D(4,3)) in the batch crystallization of sugar cane in pilot-scale process. The mathematical model presented in this work includes the population balance equation (PBE), the mass and energy balances, and the kinetics equations of nucleation and growth rate. The kinetic parameters were calculated from optimization using experimental data obtained from a pilot-scale process. An uncertainty analysis was performed and used to specify robust agitation trajectories that minimize the variations of crystal size from batch to batch. Four cases studies are presented to obtain 920, 1000, 1200, and 1300 μm of D(4,3) subject to a constraint in the formed crystal mass (FCM) of 4700 g under uncertainty in the kinetic parameters. The resulting robust agitation trajectories were implemented in the pilot-scale process. Comparisons between experimental data and the model predictions are presented. 1.3614 g/cm3). In this stage, seed crystals are introduced into the batch to induce the growth at low supersaturation levels and also to reduce the effect of spontaneous nucleation, which often produces undesirable CSD resulting in a loss of product quality.8 Batch crystallization is typically operated based on the supersaturation (Sr) of the mother liquor, which is known as the driving force for the nucleation, growth, and agglomeration of crystals. CSD and FCM can be optimized if Sr is kept into the so-called first metastable zone located between the solubility and the nucleation curve.3 Since Sr is a direct function of the crystallizer’s operating temperature, the slurry’s temperature is often used as the manipulated variable to maintain the system in the supersaturation zone. Most of applied control approaches have been performed through the implementation of time-dependent temperature profiles that aim to maintain the system in the first metastable zone. The temperature profiles are calculated based on the use of trial and error experimentation,9,10 mathematical modeling analysis,11−13 and direct design approach.14,15 Despite the efforts made in this field, the optimal operating conditions for the key variables involved in the sugar cane industrial process are not completely known and remain as an active area of research. Although most of the studies have focused on the specification of temperature and vacuum pressure profiles for the optimal operation of this process,10,16,17 the agitation rate in the system has been shown to have a significant effect on the CSD since it can potentially affect crystal growth, inhibit the spontaneous nucleation, and satisfy tight requirements in crystal size.16−21 To account for this effect, an appropriate mathematical model describing the

1. INTRODUCTION Batch crystallization is a unit operation widely used for the production of chemical, pharmaceutical, agrochemical, and food products.1−3 This operation offers several advantages over other separation techniques such as low operation costs, high purity (crystalline) products in a single (solid) stage, and attractive final product appearance for commercial purposes. The performance in batch crystallization is usually measured in terms of the product properties such as the crystal size distribution (CSD), crystal habit, and crystal morphology.4,5 CSD is of primary importance since it affects the downstream operations of separation, i.e., filtration and drying, and the product quality, which is measured by the crystal average diameter (% volume D(4,3)) and its corresponding standard deviation S(4,3). These crystal properties are directly affected by a few operational aspects such as the time at which the seed crystal is introduced into the batch, its corresponding characteristic CSD, and the cooling temperature trajectory implemented on the batch. The formed crystal mass (FCM) is another performance measure used in batch crystallization that depends on the cooling temperature profile, the vacuum pressure profile, and the properties of the initial crystal seed used in the batch. Both the FCM and the CSD are dominated by physicochemical phenomena such as nucleation, growth, and agglomeration of crystals that arise from the presence of a continuous phase and a dispersed phase.6 These phenomena, which occur at the system’s microscopic scale and that determine the quality of the product obtained from batch crystallization, cannot be accurately controlled from batch to batch and therefore lead to product variability.7 Industrial crystallization of sugar cane is divided into two sequential stages. At the first stage, the crystallizer is filled with a solution containing dissolved sucrose, also referred to as the mother liquor. This solution is then concentrated by evaporation to obtain a predefined saturation concentration (1.3594− © 2014 American Chemical Society

Received: Revised: Accepted: Published: 13180

May 1, 2014 July 9, 2014 July 26, 2014 July 27, 2014 dx.doi.org/10.1021/ie501800j | Ind. Eng. Chem. Res. 2014, 53, 13180−13194

Industrial & Engineering Chemistry Research

Article

Table 1. Devices of the Experimental Pilot-Scale Process quantity

device

2 2 1 1 6 6 1 1 2 1 1 1 1 1 1 1 1

J type thermocouple, temperature of 0−760 °C, cable length 1 m. Thermowells. Stainless steel. Vacuum pump Felisa FE-1400, 0.3 HP. Proportional control solenoid valve, Burkert. Average temperature of 140 °C, with digital controller. 2-way solenoid valve, normally closed, steel, Parker. Flow valve to allow water and steam flow through pipes. Pressure regulator Norgren Mexico. Vacuum transmitter. Cole-Parmer Model 07356-11. Stainless Steel. Pressure transmitter, Cole-Parmer Model 68072-08. Steam boiler, model MBA9 of SUSSMAN; maximum pressure, 100 Psi; work voltage, 240 VAC; control voltage, 120 VAC. Hydraulic pump, QB60 Clean Water Pump, 1.5 HP. Galvanized pipe system for water circulation. High temperature insulation system. Manometer, ASHCROFT. Condenser, stainless steel surface. Plastic tank 1100 L capacity. Stainless steel crystallizer of 12.77 L, heating−cooling jacket of 11.10 L, four vertical wall baffles of 17 cm (wide) by 3.5 cm (length), agitation arrow of 39 cm (length). Agitator for closed tank, model NSDB of HP, direct transmission of 1750 rpm (1 phase, 60 cycles), 110 VCA totally closed, without ventilation, stainless steel 316 with bridle of 4 in. (diameter) in stainless steel, with agitating arrow of 26 in. (length) and 1/2 in. of diameter in stainless steel 316; velocity investor (driver) integrated with rank from 0 to 1750 rpm Marine propeller type impeller. Programmable tachometer. Range from 50 to 999,990 rpm. Optical sensor for distances of 3 ft. Range 1−150000 rpm. Microscope trinocular 48923−30 Cole-Parmer. Monochrome camera with RS-170, video lens with 0.19 mm per pixel. PC, Intel Pentium IV. Operating system XP, 4 GB RAM memory, hard disc of 1 TB. Data acquisition card (NI PCI-6023E, NI PCI-6025E, NI PCI-6711, and NI PCI-232/2). Analogic−digital and digital−analogic converters allow the input/ output analogic and digital signals. Signal conditioning module (NI SCC-TC02). Image acquisition card (NI PCI-1409). Galvanized steel condenser of direct contact with 2.27 m high and 0.3 m of diameter.

1 1 1 1 1 1 1 4 2 1 1

(FCM) or poor crystallization (excessive nucleation or very small crystals). Accounting for uncertainty while designing optimal temperature profile and agitation rate trajectories is therefore essential for maintaining the process within feasible operating regions where crystallization can be successful. Batch crystallization is often modeled using the population balance equation (PBE), which allows the derivation of a system of nonlinear partial differential equations that describe the accumulation rate of the CSD.29,30 The PBE describes the evolution of CSD, from which the moments of distribution can be deduced and can be coupled with the system’s mass and energy balances to describe the accumulation rate of the state variables, e.g., temperature, vacuum pressure, and concentration. The complex nature of the mathematical models of particulate processes has motivated extensive research efforts on the development of numerical methods and/or accurate analytical solutions for the PBE.11,31−33 The method of moments (MOM)26,27 is currently one of the most used mathematical techniques to solve the PBE, where only loworder moments of the CSD are calculated and unknown parameters of an assumed distribution are adjusted to the calculated moments. A limitation of the MOM is that the procedure must form moments from the PBE, which may lead to cases where the resulting equations may require considerable analytical effort or cannot be treated with this method. This work presents a study on the effect of agitation rate on the CSD and FCM in a vacuum batch crystallization of sugar

time evolution of the CSD with adequate nucleation and growth kinetics9,22 and that accounts for the agitation rate effects is needed. Although studies that included the agitation rate as a manipulated variable have predicted improvement on the CSD, the optimal trajectories calculated by those analyses have shown a decrease in effectiveness in real applications.23,24 This loss in performance may be associated with the uncertainty in the kinetic parameters of model, which arises in (a) microscale and (b) macroscale. The microscale events occur to a lesser extent, and they are associated with stochastic events in the system such as the crystal’s nucleation and growth. The macroscale events are associated with time-varying parametric uncertainty, i.e., physical parameters that are unknown and may slowly change their values during the evolution of the process,25−27 external perturbations, and, e.g., impurities in the feedstock solution. This work considers as a key innovation that uncertainty in the kinetic parameters is implicitly assumed into a time-varying deterministic model and sensitivity analysis rather than only being associated with microscale events. The uncertainty associated with the macroscale events affect the CSD, and it depends on the width of the concentration zones (first and second metastable zone) and supersaturation value of the process at a given time. The width of the concentration zones may vary significantly (density > 1 × 10−3) depending on the operating temperature28 and the agitation rate.24 The latter may result in either satisfactory crystallization: large crystal size (D(4,3)) and high amount of crystal mass 13181

dx.doi.org/10.1021/ie501800j | Ind. Eng. Chem. Res. 2014, 53, 13180−13194


Article

Figure 1. Experimental batch crystallizer (pilot-scale process).

2. EXPERIMENTAL SET-UP

cane into a pilot-scale process. The key novelty introduced in this work is that uncertainty in the kinetic parameters is implicitly assumed into a time-varying deterministic model, the optimal design, testing, and implementation of robust agitation rate trajectories on the pilot-scale process. The agitation rate trajectories were obtained from optimization using a PBE model developed for a pilot-scale operation of sugar cane. The uncertainty in the kinetic parameters was obtained from experimental data generated from an actual pilot-scale process. The results obtained from the implementation of the optimal agitation rate trajectories on the pilot-scale process are presented and discussed in this study. To the authors’ knowledge, this is the first study that computes and implements optimal agitation rate trajectories obtained from a PBE model that explicitly considers uncertainty in the system’s kinetic parameters. This article is structured as follows: the next section presents the experimental setup (pilot-scale process) used in the experimental stage to obtain the necessary experimental data for the kinetic parameter estimation and experimental validation. The mathematical framework used in this work and the novel strategies used to calculate the rate trajectories of optimal agitation while taking into account the uncertainty in the process is presented in the Methodology section. The Results and Discussion section presents the kinetic parameters calculation, uncertainty analysis, optimization of agitation rate trajectories, case studies of crystal size D(4,3), and finally experimental validations and remarks. Conclusions are presented at the end of this article.

The equipment used in the experimental stage (detailed in Table 1) is a stainless steel batch crystallizer (pilot-scale process) with a heating−cooling jacket, steam generator, DC motor, vacuum pump, a direct contact condenser, and a tachometer that specifies the agitation rate (see Figure 1). The study reported by Bolaños et al.10 provides a complete description of all the equipment and instrumentation devices mounted on the pilot-scale process of sugar cane. The operation of this batch crystallization pilot-scale process is described next. The operation begins with the loading of the saturated sugar solution (1.3594−1.3614 g/cm3 at 70 °C) into the crystallizer. In the case the solution is unsaturated, solvent (water) must be evaporated to reach the saturation concentration by adding steam into the heating−cooling jacket. Vacuum pressure is set to 76.20 kPa and is kept constant during 40 min; the agitation rate is usually set to 300 rpm. In the first minute of the process, the crystal seed is introduced into the crystallizer. After 40 min (76.2 kPa of vacuum pressure), the vacuum pressure is increased until 84.66 kPa following a predefined natural cooling temperature profile. Sampling is performed in the process to obtain information from the experimental concentration (slurry) and CSD at each sampling time (15 min). The methodology to perform this activity is as follows: 1. The sampler device (SD) is introduced vertically into the batch (see Figure 1) to obtain 10 mL of slurry (crystals 13182



Article

rate, and the parameters D(L) and B(L) represent the death rate and birth rate, respectively, which arise due to the agglomeration of the crystals or breakage; these terms are often represented as α(L), which is also known as the production− reduction term. To simplify the mathematical model, the following assumptions are often considered for the PBE model:36,37 i. Agglomeration and breakage of crystals is negligible. ii. The secondary nucleation rate is considered as the total nucleation rate. iii. The crystal nuclei produced have negligible size. iv. The system is well mixed. v. The crystals inside the tank are well suspended, i.e., no accumulation of crystals at the bottom of the crystallizer. Based on the above, the MOM is applied to eq 1 producing the following set of ordinary differential equations (ODEs):36,37

and solution) without breaking the equipment’s vacuum pressure. 2. A separation is performed on the sample through filtration to obtain a solution without crystals. The density of the sample is measured in the DMA 4500 device (dissolved mass crystals per unit volume solvent) (see Table 1). The crystals obtained from the filtration process are used to measure CSD by using imaging analysis. The procedure to quantify CSD is described next. 2.1. Image Acquisition System and CSD Calculation. IMAQ Vision Builder (National Instruments, Inc.) was used to analyze the crystals obtained at every sampling time of the experimental runs, i.e., every 15 min. This image-based approach is an alternative in measuring both length and area of particles (crystals) in a direct way through IMAQ Vision’s software. This technique consists of acquiring an image using a monochromatic camera with video RS-170 and 60 Hz crisscross (8 bits of resolution) and handling the light beam from a microscope trinocular. The sample is transported from the crystallizer to the imaging acquisition system by means of a peristaltic pump running at 1500 rpm. The camera captures 10 square images (which includes approximately 100 crystals) that are processed and cleaned, which avoids undesirable light variations and is achieved through the threshold technique that allows the specification of a single image in the gray scale. Areas with high density of crystals are isolated to be independently analyzed, and black pixels (crystals) are counted. The black pixels are compared with specific standards to decide if the crystal is present or not according to binary images (background).34 A Neubauer’s camera was used to count the particles and determine a conversion factor through simple calibration, i.e., the conversation factor enables a direct relation from 1 pixel (one pixel side) to 200 μm (length). In this work, a microscope with a 10× ocular lens, a 40× objective and an E square from Neubauer’s camera from 50 μm away was used for the measurements and analysis of the particles. This is equivalent to 20 000 μm (10 × 40 × 50) per 100 (length of a pixel side).10 Thus, 1 μm is equal to 0.005 of the length of a pixel side.

dμ0 dt

dμk dt

k = 1...5

(4)

B0 = KbSrbMTj Nrp

(5)

G = K gSrg Nrq

(6)

Sr =

C − Csat Csat

(7)

The aim of this study is to determine accurate uncertain descriptions for the kinetic’s power-law-type uncertain parameters (i.e., Kb, b, j, p, Kg, g, q), optimizable from experimental data. Thus, this work does not consider the mechanistic approach (i.e., frequency and energy of collisions between particles with the impeller) of the crystallization kinetics. The application of the kinetic approach using powerlaw-type has recently been used to study different processes of crystallization at medium- and large-scale levels.38−41 From eq 5, MT represents the total mass: M T = ρc K vμ3 (t )

(1)

(8)

where the third moment μ3(t) is represented as

n(L , t ) = n0(L , t )

μ3 (t ) =

t=0

∫0

∞

nL3 dL

(9)

The system’s mass balance is

and the boundary condition B0 GL = L0

= kGμk − 1(1 + ln(V ))

where Sr represents the relative supersaturation defined by

with the initial condition

n(L0 , t ) =

(3)

Only the first four moments of the distribution have a physical meaning, i.e., number of particles, length, area, and volume. The nucleation (B0) and growth (G) rates are modeled with empirical power law type equations defined by9,22

3. METHODOLOGY 3.1. Mathematical Framework. The PBE considered in this work assumes that growth is independent of the crystal size35 and that changes in volume (continues phase) are due to vacuum evaporation effects. The PBE model is as follows: V ∂(Gn) ∂(nV ) = V (α(L)) + ∂L ∂t

= B0

d(C) = −ρc K vh[3 dt

(2)

∫0

∞

GnL2 dL + B0 L0 3]

(10)


where the term ∂n/∂t represents the change in the number of density with respect to time in a batch crystallizer, ∂(Gn)/∂L describes the difference between crystal growing inside and that growing outside a range of crystal size dL due to crystal growth

C(0) = C0

(11)

Similarly, the energy balance inside the crystallizer is16,39 13183


Industrial & Engineering Chemistry Research ρVCp

Article

dT dPV = c − HevFev − ΔHcρc dt dt K vV [3

∫0

∞

profile while considering different agitation rate trajectories. Thus, the present analysis assumed that both p and q were constant parameters and equal to a nominal value. On the other hand, those studies also observed that b and g (for Sr) in eqs 5 and 6 have the main higher effect on the batch crystallization process for different cooling temperature profiles. The selection of a natural cooling temperature profile over other profiles, e.g., cubic and linear, was based on the results presented by Bolaños et al.10 that showed improvement on the growth of D(4,3) and FCM when using that temperature profile. Based on the above, the following least-squares optimization problem was formulated where the sum of squared errors (SSE) between the experimental data and the crystallizer model (MOM model) is minimized for each sampling time i in the batch, and then

GnL2 dL + BoL0 3]

+ U1A1(Tj − T )

(12)


T (0) = T0

(13)

where Kv = π/6 is the characteristic shape factor of the sugar crystal,42 whereas ΔHc is the crystallization heat, i.e.,43 ΔHc = −12.2115 − 0.7937T

(14)

Dynamics of liquid volume (continues phase) inside the crystallizer due to the evaporation of solvent (water) are described with the expression shown in eq 15. That expression is considered for a vacuum pressure set to 76.2 kPa, with 40 min of total evaporation time (constant). Then the process follows a vacuum pressure profile (type natural cooling) in order to force by resetting of crystallizer’s vacuum pressure to 84.66 kPa.44 Table 2 shows the complete set of parameters used in the present analysis. V = 752.3333 − 24.271t + 0.3098t 2 − 0.00145t 3

⎛ CSDexp, i − CSDmodel, i ⎞2 ⎟⎟ min ∑ ⎜⎜ Ψŝ CSD ⎝ ⎠ i exp, i=1 n

subject to : Crystallizer model eqs 3−15

where Ψ̂s represents the nominal vector of the kinetic model parameters estimates, i.e., Ψ̂s = [Kb, b, j, p, Kg, g, q]. CSDexp and CSDmodel represent arrays with D(4,3) and FCM data obtained from the actual pilot-scale process (Figure 1) and the MOM model at particular times, respectively. Note that the errors have been normalized in the problem’s objective function to equally weight the deviations in D(4,3) and FCM, respectively. This was considered to obtain nominal estimates of the kinetic parameters that can predict large crystals in a large volume. Problem 16 was solved using the constrained nonlinear optimization programming function ( f mincon) available in MATLAB. The set of ODEs representing the crystallizer’s mathematical model was solved using the ode23s built-in function, which implements a modified version of the second order Rosenbrock implicit method.45 3.3. Uncertainty Analysis. In batch crystallization, reproducibility is a frequent problem4 that has significant effects on the quality of batch-to-batch runs and that is mostly caused by time-varying parametric uncertainty, i.e., physical parameters that are unknown and may slowly change their values during the evolution of the process,25−27 external perturbations, e.g., impurities in the feedstock solution, and to a lesser extent, the microscale (stochastic) events occurring in the system such as the crystals’ nucleation and growth.46−48 In order to account for product variability, the present analysis performed first a sensitivity analysis to identify the kinetic parameters of the model that have a significant effect on the crystal average diameter (% volume D(4,3)), which is the most important product property in batch crystallization of sugar cane, and the FCM, which measures the economic performance of this process. The sensitivity analysis was performed through the variation in ±5, 10, and 20% of the nominal parameter estimates obtained for Ψ̂s from the solution of problem 16. The sensitivity analysis showed that the nucleation rate’s parameters Kb and p and the growth rate’s parameters Kg, g, and q are the most sensitive parameters to D(4,3). However, these parameters do not result in sensitivity for the FCM (where the observed variability was much less than observed variability to D(4,3)). To simplify the analysis, these kinetic parameters were lumped in a single vector, i.e.,Ψs = ⌊Kb, p, Kg, g, q⌋. An analysis aimed at determining the amount of variability that needs to be included on the system’s most sensitive kinetic

(15)

Table 2. Process Parameters for Solving Mass, Energy, and Population Balances variable

value

units

ρc Kv A1 L0 Vj U1 ρw Cpw Cps T0

1.588 π/6 2004.4648 0.0019126 11102.9 0.3943 8.02 0.5862 2.4687 70

g/cm3 -cm2 cm cm2 J/(°C·min·cm2) g/cm3 J/g·°C J/g·°C °C

(16)

3.2. Computation of the System’s Kinetic Parameters. As shown in eqs 5 and 6, the proposed mathematical model includes seven kinetic parameters, i.e., Kb, b, j, p, Kg, g, and q, which are specific for each operating condition in batch crystallization. The values of these parameters need to be adjusted to ensure that the mathematical model provides representative results of the actual pilot-plant process. In this work, nominal estimates for the kinetic model parameters were obtained from experimental data from the batch crystallization system presented in Section 2. These experimental data were previously reported by Alvarado44 and were obtained under the following operating conditions: • Total evaporation time of 40 min with 76.2 kPa of vacuum pressure (58 °C). • Natural cooling temperature profile (adiabatic cooling (58 to 41 °C) forcing a change in the vacuum pressure from 76.2 to 84.66 kPa). • Batch time of 90 min. • Agitation rate of 225 rpm According to Bolaños16 and Quintana et al.,9 the parameters p and q (for Nr) in eqs 5 and 6 do not change significantly (i.e., by an order of magnitude) for the natural cooling temperature 13184



Article

work that has not been considered before for batch crystallization. 3.4. Optimal (Robust) Agitation Rate Trajectories under Uncertainty. Agitation rate plays an important role in batch crystallization.16,21 At the final batch times, low agitation rates produce temperature and concentration gradients inside the system that may lead to the formation of local nucleation and growth rates, which result in different CSD, thus affecting the process performance. Quintana et al.9 suggested that the natural cooling profile, in combination with linear agitation rate trajectory, produced the best estimates for D(4,3) (507 μm) and FCM (626 g). Likewise, Alvarado44 showed that a natural temperature cooling profile at a constant agitation rate of 225 rpm produced a D(4,3) of 1149.53 μm and FCM of 5003.0 g, which is a satisfactory D(4,3) since large crystals are produced that can be completely recovered using centrifugation in the downstream separation processes, thus preventing the loss of product. Based on the above, this study considers the natural cooling temperature profile as one of the most suitable manipulated variables for the crystallization process. On the other hand, the implementation of a suitable agitation rate profile may further improve the process dynamics and therefore lead to improvements in product quality and process economics. In addition, product variability can also be explicitly accounted for by incorporating the uncertainty effects of the system’s kinetics in the analysis. This results in the specification of robust agitation rate trajectories that can achieve particular (user-defined) product quality properties in the presence of uncertainty in the system’s kinetic parameters (Ψs). Based on the above, optimal (robust) agitation rate trajectories can be obtained from the next optimization problem:

model parameters is next performed using the observed variability in D(4,3) from experimental data. The experimental data from FCM was considered in optimization problem 17 to ensure a maximum limit of produced crystal mass, due to kinetic parameters of the model (Ψ̂s) that are function of a finite crystal mass (by mass balance), since this process was operated in batch mode. n

min

δΨ∈ s [0,1]

∑ (eu ,i + el ,i) i=1

subject to : Crystallizer model eqs 3−15 eu , i

⎛ D(4, 3)model, u , i − D(4, 3)MaxExp, u , i ⎞2 ⎟⎟ , = ⎜⎜ D(4, 3)MaxExp, u , i ⎝ ⎠

⎛ D(4, 3)model, l , i − D(4, 3)MinExp, l , i ⎞2 ⎟⎟ , el , i = ⎜⎜ D(4, 3)MinExp, l , i ⎝ ⎠

i = 1...n

i = 1...n

Ψs = Ψŝ ± Ψŝ δ Ψs (D(4, 3)model, u , i − D(4, 3)MaxExp, u , i ) ≥ 0, (D(4, 3)model, l , i − D(4, 3)MinExp, l , i ) ≤ 0, FCM model, u , i ≥ FCMMaxExp, u , i,

i = 1...n i = 1...n

i = 1...n (17)

Solution of problem 17 gives the region of experimental variability (δΨs) for the kinetic parameters around their nominal values (Ψ̂s), at each sampling time (i = 1, ..., n); where eu,i and el,i represent the normalized (squared) deviations in the positive (maximum) and in the negative (minimum) direction between the model predictions and the experimental data reported for D(4,3). D(4,3)MaxExp,u,i and D(4,3)MinExp,l,i represent the maximum and minimum experimental values observed for D(4,3) at each sampling time i. The FCMmodel,u,i and FCMMaxExp,u,i represent the maximum predicted values from the crystallizer model and maximum experimental value at each sampling time i. Note that multiple experimental data sets were used in the present analysis to account for the system’s dynamic variability. The optimization problem shown in eq 17 proceeds as follows: for each set of values in δΨs selected by the optimization algorithm, the crystallizer model equations are solved for each combination between the extreme values of the uncertain kinetic parameters included in the vector Ψs. This results in a set of D(4,3) estimates for each sampling interval i. Then, the minimum and maximum expected D(4,3) values at each time interval, which are due to the combination between the extreme limits of the uncertain kinetic parameters, are identified (i.e., D(4,3) model,l,i and D(4,3)model,u,i). These estimates are then compared against the maximum and minimum D(4,3) observed from the experimental data, i.e., D(4,3)MaxExp,u,i and D(4,3)MinExp,u,i. The result of this comparison is used to evaluate the problem’s objective function. The optimization problem 17 returns bounds on the most sensitive kinetic parameters that affect the production of crystals obtained from sugar cane batch crystallization. Note that the bounds obtained for the uncertain kinetic parameters have been calculated from experiments. This is a key novel aspect of this

min

rpm ∈ [50,700], t f

tf

subject to : Crystallizer model eqs 3−15 D(4, 3)l |t = t f ≥ D(4, 3)target FCMl|t = t f ≥ FCM target t = [0, t f ]

(18)

where D(4,3)l|t=tf and FCMl|t=tf represent the lower predicted values from the crystallizer model and D(4,3)target and FCMtarget represent user-defined product targets that need to be satisfied at the final batch time tf. The optimization problem 18 was solved following the same methodology of problem 17 for the selection of δΨs. This results in a value of D(4,3)l|t=tf estimated to be compared with D(4,3)target to satisfy the minimum D(4,3) specified for quality propose. Note that eqs 3−15 describing the transient behavior of the pilot-scale process of sugar cane batch crystallization explicitly account for the uncertainty in the kinetic model parameters (Ψs) identified from problem 17. As shown in problem 18, an optimal (robust) agitation rate trajectory ensures that the specifications for D(4,3) and FCM are satisfied at the final batch time in the presence of uncertainty in the system’s kinetic parameters. Although another product’s properties can be added into this problem (e.g., S(4,3)), D(4,3) and FCM are by far the most 13185



Article

Table 3. Kinetic Parameters nucleation kb, No. Part./cm3·min·(g/cm3)j·(rpm)p −2

1.56 × 10

growth

b

j −4

6.27 × 10

p −3

2.04 × 10

representative properties used to measure the quality of sugar crystals at the industrial scale. Hence, problem 18 aims to specify the agitation rate profile and the smallest final batch time needed to achieve specific properties of the crystals obtained from this process in the presence of uncertainty in the system’s kinetic parameters. Problem 18 was solved for four scenarios that involve different D(4,3) targets (i.e., 920 μm, 1000 μm, 1200 μm, 1300 μm) and a FCM production of at least 4700 g at the final batch time. To study the agitation effects on the D(4,3), the batch time was discretized at intervals of 20 min in which agitation rate can vary in a range between 50 to 700 rpm. In a large-scale crystallizer, these changes in rpm are possible because these are only mechanical changes and not phenomenological. Also, these changes involve a lag time about 3−10 s, which is not significant in large-scale sugar cane batch crystallization.49 In contrast to the methodology proposed by Chianese and Kramer21 where the problem’s cost function is based on the D(4,3) and batch time remains constant during the optimization, the present work aims to search for the robust agitation rate profile that achieves a desired D(4,3) while reducing the pilot-scale process’ operating costs, i.e., at shorter batch times. The present work is the first study that explicitly accounts for uncertainty in the system’s kinetic parameters calculated from experiments, i.e., the amount of uncertainty considered in the specification of the robust agitation rate profiles from problem 18 was not arbitrarily specified but has been estimated from experimental data. These aspects are key novelties introduced by the present work that have not been previously reported in the open literature. Problem 18 can be casted as a nonlinear dynamic optimization problem. To ensure that the trajectories obtained from problem 18 are near optimal, this work made use of the function multistart in MATLAB, which performs parallel optimization varying the initial value (randomly) of the decision variables. An analysis was performed to determine the necessary number of iterations needed to ensure the reliability of the results without increasing the computational time. It was concluded that 40−50 iterations are appropriate since no reduction in the cost function was observed when using a larger number of iterations. Furthermore, the scenarios considered in this study were obtained at the operating conditions described in Section 3.2.

1.41

kg, cm/min·(rpm)q 1.31 × 10

−05

g

q

1.035

8.45 × 10−1

Figure 2. MOM model results: (a) D(4,3), (b) FCM, (c) slurry temperature, (d) nucleation rate, and (e) growth rate.

from problem 17 and the experimental data reported by Alvarado44 are in reasonable agreement; i.e., the MOM model predicts the experimental data with maximum deviations of 50 μm. At time 40 min, which is the time at which the natural temperature cooling profile begins, it can be observed that the MOM model is sensitive to changes in the vacuum pressure by increasing the value of D(4,3). Figure 2b shows the FCM behavior; in this case, prediction is adequate obtaining approximately the same FCM from both the MOM model and the experimental data at the end of the batch. Figure 2c shows that the MOM model predicts well the slurry’s temperature, which is presented in three stages. Stage 1 comprises the first 5 min in the batch crystallizer operation, where the system reaches the vacuum pressure for the constant evaporation stage (stage 2) and where the MOM model represents the sudden cooling properly to reach a temperature of 58 °C; vacuum pressure then remains constant at 76.2 kPa during 40 min, where large amounts of solvent (water) are evaporated. Further, the natural cooling profile is implemented by introducing a step change in the system’s vacuum pressure at stage 3. As shown in Figure 2c, the MOM model accurately predicts the cooling in the system reaching a final temperature of 41 °C.

4. RESULTS AND DISCUSSION 4.1. Kinetic Parameters. Table 3 shows the nominal estimates for the system’s kinetic parameters obtained from problem 17. These estimates were calculated for a natural cooling profile, a total evaporation time of 40 min, and vacuum pressure of 76.2 kPa. It is important to establish that the obtained nominal estimates for the nucleation parameters are small because the crystallization process considered was seeded to reduce the nucleation of new crystals; this also promote that process starts above the saturation line (into first metastable zone) where the seed did not have dissolution. As shown in Figure 2a, the calculated fitting for D(4,3) from the MOM model using the nominal kinetic parameters’ estimates obtained 13186



Article

Table 4. Uncertainty in Kinetic Parameters nucleation rate δkb 48.301 × 10

growth rate δp

−3

δkg −3

72.743 × 10

δp −02

12.054 × 10

δq −02

18.327 × 10

28.799 × 10−3

Figure 3. D(4,3) variability, FCM variability.

literature.19−21,24 Figure 3 shows the experimental D(4,3) and the expected D(4,3) variability calculated using multiple combinations between upper and lower bounds of kinetic parameters. As shown in Figure 3, the region delimited by the D(4,3)’s variability includes the experimental data. These results allow possible range of values for D(4,3) during (and at the end) of a batch to be predicted, which enables the standardization of product quality. The experimental data points of D(4,3) lay in the shaded region (D(4,3) variability) because they were used in the formulation of problem 17 and are expected to fit that region. On the other hand, the shaded region of FCM variability was obtained as a prediction result from the crystallizer model (3−15) and problem 17. This only represents an approximation of FCM experimental data as shown in Figure 3. According to the formulation shown in problems 17 and 18, the uncertainty specified for each kinetic parameter is quite small (Table 4), i.e., the uncertainty represents 0.048%, 0.072%, 0.12%, 0.18%, and 0.028%, for Kb, p, Kg, g, and q, respectively. This novel strategy has never been applied before for batch crystallization and is a key contribution of this paper. The five uncertainty parameters (Kb, p, Kg, g, and q) are small in variability over its nominal value (Table 3) because they have high sensitivity; thus, changes in those parameters by a small amount (Table 4) produces high variability in the key process variable, i.e., the crystal’s size D(4,3). 4.3. Optimal Agitation Rate Trajectories. 4.3.1. Case Study 1: D(4,3)target = 920 μm. According to Alvarado,44 the crystal obtained at the end of batch run while using the operating conditions reported in that study has a D(4,3) of approximately 1000 μm. Under those conditions, the minimum D(4,3) that the MOM models predicts (using the nominal estimates for the model’s kinetic parameters reported in Table 3) is 907 μm with a constant agitation rate of 50 rpm. Therefore, the aim of the present case study is to search for an optimal agitation rate trajectory that achieves a D(4,3) of 920 μm in the shortest possible batch time. The agitation rate profile was initially obtained using the nominal kinetic parameter estimates depicted in Table 3. Therefore, problem 18 was solved under perfect knowledge of the system’s kinetic

Figure 2a,b also shows that there exists inherently variability in the experimental data from batch-to-batch runs although they have been made under the same operating conditions. This is a clear indication of the need to account for uncertainty in the crystallization process. Figure 2d,e shows the estimated nucleation and growth rates from eqs 1 and 2 obtained while using the kinetic model parameters shown in Table 3, respectively. Due to the difficulty in measuring these rates, experimental data are not available for these kinetic variables. For the growth rate, there are two characteristic stages; the first stage is between 0 and 40 min (constant evaporation) where growth is kept constant. After this stage, the natural cooling of the solution promotes an initial increase in crystal growth and subsequently decreases toward the end of the batch. For both kinetics (nucleation and growth), the behavior is similar to that observed by Quintana et al.9 and Bolaños.16 4.2. Uncertainty Analysis. The results shown in Figure 2 suggest that there exists variability in the system and that the discrepancies between the mathematical model used here to describe the batch crystallization pilot-scale process and the actual pilot-scale process data can be better captured by accounting for uncertainty in the system’s kinetic parameters. As discussed in Section 2, a sensitivity analysis was performed on the nucleation and growth rate kinetic parameters to identify the most sensitive kinetic parameters, i.e., Kb and p from the nucleation rate and Kg, g, and q for the growth rate. Table 4 shows the amount of uncertainty obtained from problem 17 that need to be considered for each of these uncertain kinetic parameters. Although the uncertainties for parameters p and q are relatively small when compared to the other uncertain parameters, these two parameters have high sensitivity to D(4,3) so that changing these parameters’ estimates by a small amount produces high variability in the crystal’s D(4,3). It is also noted that the effect of the uncertain kinetic parameters on D(4,3) is different, i.e., increasing the nominal value of kb, p, and g produces a decrease in growth toward D(4,3) whereas increasing kg and q promotes growth in D(4,3). Those effects were tested using the mathematical framework developed in this work, and the results agree with those reported in the 13187



Article

Figure 4. Nominal agitation rate profiles. (a) case study 1, (b) case study 2, (c) case study 3, and (d) case study 4.

time. The results obtained for the present case study agree with those reported in the literature.16,20 The calculation of the optimal agitation rate profiles (both constant and time-varying) while using uncertain kinetic model parameters is next considered. Figure 5a shows the results obtained from that analysis. The batch times needed to achieve the desired product quality target were 95.35 and 115.7 min for the cases of time-varying and constant agitation rate profiles, respectively. These results indicate a reduction in batch time of approximately 17.55% in favor of the time-varying robust agitation profile. By comparing Figures 4a and 5a, it can be observed that the optimal agitation rate trajectories obtained with and without uncertainty in the system’s kinetic parameters are relatively different. In the interval of 40−60 min the step decrement in the robust agitation rate profile is performed in a moderate way in comparison with the nominal agitation rate case, which results in larger batch times. For both nominal and robust cases, the agitation rate followed a somewhat similar

model parameters and equal to their nominal estimates listed in Table 3. To compare this result, problem 18 was also solved for the case of constant agitation rate profile. In that case, the goal was to search for the optimal (constant) agitation rate that satisfies the product quality constraints while using the shortest batch time. Figure 4a shows the results for D(4,3), FCM, and the agitation rate trajectories obtained while using nominal estimates for the kinetic parameters. As shown in that figure, a nominal agitation rate trajectory with a batch time of 90.68 min and a constant agitation rate (60 rpm) with a batch time of 117 min were obtained from this case study. That is, a reduction of 22.48% in batch time is achieved when using a nominal timevarying agitation rate profile. As shown in Figure 4a, the agitation rate is initially set high to promote faster nucleation and, therefore, creation of crystals; then, the agitation rate is suddenly reduced at time = 40 min to promote the growth of the crystals produced in the first stage of the batch until a desired D(4,3) and FCM are reached at the end of the batch 13188



Article

Figure 5. Robust agitation rate trajectories: (a) case study 1, (b) case study 2, (c) case study 3, and (d) case study 4.

path in D(4,3) and FCM, respectively. These results enable the identification of a variable effect of the agitation rate over the D(4,3), i.e., for D(4,3) with less than 1000 μm, it is desirable to start the batch with an intense agitation rate (600−700 rpm) whereas for D(4,3) greater than 1000 μm, the agitation rate should be operated between 50 and 100 rpm (see the discussion in Section 4.3.3). This effect is similar to that reported in Chianese and Kramer,21 where the effect of the agitation rate is variable with respect to the batch time. 4.3.2. Case Study 2: D(4,3)target = 1000 μm. This case study aims to specify optimal agitation rate profiles that satisfy a D(4,3) of 1000 μm and a FCM of 4700 g at the end of the batch. As shown in Figure 4b, the batch time required for the case of time-varying and constant agitation rate (604 rpm) profiles is the same while using the nominal kinetic parameter estimates in the solution of problem 18, i.e., approximately 79

min. The same behavior was observed when uncertainty in the kinetic model parameters was considered in the computation of the time-varying and constant agitation rate profiles, i.e., the batch time reported for both cases is roughly 85 min. Note that in the case of the constant agitation rate profile, the agitation rate needed to increase the crystal’s D(4,3) to 1000 μm is almost seven times the rate needed to produce crystals with a D(4,3) of 920 μm (see Case Study 1). The application of high agitation rate throughout the batch for the constant case (667 rpm) can generate crystal breakage, decreasing the performance in real applications and the lifetime of the agitation system. 4.3.3. Case Study 3: D(4,3)target = 1200 μm. This case study was conducted to examine the possibility of increasing the D(4,3) to 1200 μm at the final batch time while producing a mass of crystals above 4700 g. Figure 4c shows the D(4,3) and FCM obtained for this case study (nominal kinetic parameters) 13189



Article

Figure 6. Sugar cane crystals from experimental validation. Representative images from case study 3.

that may generate crystal breaking, thus decreasing the system’s performance (desired D(4,3)). The reduction in batch time is 25.58% when compared to the constant agitation rate profile. A particular case occurred in the robust case: a constant agitation rate that satisfies the requirements of D(4,3) and FCM could not be obtained from problem 18. For that reason, Figure 5d shows results only for the robust (time-varying) agitation trajectory. This result shows that specifying optimal timevarying agitation rate profiles enables the production of high purity crystals with larger crystal sizes, a property particularly desired for batch crystallization products. High agitation rates (600−700 rpm) were used in the case studies for the optimized agitation rate trajectories for nominal and robust cases (Figures 4a,b and 5a,b); these may have a negative effect in the process, e.g., crystal breaking, agglomeration, and dissolution. However, these aspects are not significant since (a) the time is no more than 20 min longer in relation to the final batch time (90−120 min), (b) high agitation rates at the beginning of the process avoid segregation in the crystallizer, and (c) no impact is observed in substances with high viscosity such as the one used in this work (sugar crystals). 4.4. Experimental Validation. At Figure 7, the results of experimental validation (run and its replicate) against the simulations (MOM model) presents a satisfactory fit (90%), which indicates that the optimization of the seven parameters involved in eqs 5 and 6 with the power-law-type model is adequate as the reproducibility of the phenomenon described is achieved. Experimental runs were conducted in the vacuum batch crystallizer (see Figure 1) to validate the nominal and the robust time-varying agitation rate trajectories obtained in case studies 1 and 3. These case studies were selected because their trajectories have an inverse effect on D(4,3) and the experimental validation allows us to conclude properly on the effects of the agitation rate. Figure 7a,c shows the results for D(4,3), FCM, and the concentration−temperature diagram when the time-varying agitation rate trajectories are obtained for the nominal case and those shown in Figure 4a,c, respectively. According to Figure 4a for D(4,3), it is expected that crystal growth is high during the first 40 min of batch time, and then the growth trend is reduced to reach 920 μm. As shown in

as well as their corresponding agitation rate profiles. In this case, the D(4,3) has a more conservative beginning when compared to the profiles obtained from the previous cases studies. This is because of the agitation rate trajectory set to the process, which begins with a slow agitation of 50 rpm; 60 min later, a step change is performed to reach 500 rpm, which leads to massive migration of the solute available in the solution to the crystal faces already formed, generating additional growth of crystals. This effect may be related to the high viscosity of the solution (sugar crystals) which reduces the mass transfer and likely the collision frequency of particles with the impeller (Figure 6). Operating times are 100.7 min for a constant agitation rate (171.89 rpm) and 86.3 min for the time-varying trajectory, resulting in a 14.24% reduction in batch time for the latter. The agitation rate profile shown in Figure 4c indicates that there is a positive step in the range of 20−40 min. This occurs due to the constraint for FCM of 4700 g, and the subsequent step changes ensure that the desired D(4,3) is achieved (1200 μm). Another consideration in the application of the agitation rate is the reduction of the energy consumption, since applying high agitation rates at the final batch time improves heat flow from the slurry to the cooling fluid in the jacket, allowing a significant decrease in the amount cooling water required. The robust case (Figure 5c) gives similar results for D(4,3) and FCM as in the nominal case, resulting in a batch time of 124.3 min for the constant agitation rate and 100 min for the time-varying trajectory, which means a reduction of 19.54% in batch time. The robust agitation trajectory is similar to the nominal case, with the difference that the step that ensures compliance with the FCM constraint is performed in the time range of 0−20 min. According to Sander and Kardum,20 the step change is expected since an increase in the agitation rate in the initial stages of the batch (0−40 min) promotes the growth of the D(4,3). 4.3.4. Case Study 4: D(4,3)target = 1300 μm. This case study was performed to analyze the maximum D(4,3) that can be achieved in the current pilot-scale sugar cane batch crystallizer with the operating conditions considered for this case study. The results are similar to those obtained for case study 3 (see Figure 4d); the difference is the agitation rate trajectory, which is aggressive in the step change performed at time 60 min and 13190



Article

Figure 7. Experimental validation. (a) Nominal agitation rate trajectory for 920 μm, (b) robust agitation rate trajectory for 920 μm, (c) nominal agitation rate trajectory for 1200 μm, (d) constant agitation rate for 1200 μm, and (e) robust agitation rate trajectory for 1200 μm.

Figure 7a, the crystal growth is intense during the first 40 min of batch time. Then, the MOM model overestimates the crystal’s D(4,3) when compared to the actual crystal’s growth, i.e., the crystal’s growth is slowed after 40 min, reaching a D(4,3) of 758.09 μm at the end of batch. This can be attributed to the existence of natural sedimentation (gravity effect) at the

bottom of the crystallizer due to slow agitation (50−100 rpm), and minor importance is given to mechanical limitation of the sampling device used in the crystallizer (vertical location), which only allows (small) crystals suspended in the solution to be measured and not the largest crystals located at the bottom of the crystallizer. The same behavior can be observed in Figure 13191



Article

7b (robust case); i.e., the experimental D(4,3) remained under the uncertain region after 60 min resulting in a final D(4,3) of 790.12 μm. In Figures 4c and 5c (case study 3), it can be observed that the D(4,3) is affected by the changes made in the agitation rate; the experimental dynamics of D(4,3) (see Figure 5b) during the batch time is similar to that predicted by the MOM model; however, there are differences when the agitation rate is low (50−100 rpm), since the model has a smaller crystal growth (200 μm) when compared to the experimental data. This phenomenon occurs because the kinetic parameters used in the MOM model were calculated from experimental data reported by Alvarado44 using a constant agitation rate of 225 rpm. For this reason, low or high agitation rates, e.g., 50 and 600 rpm, can present deviations from the experimental data. According to Akrap et al.50 and Ni and Liao,24 the concentration zone widths are modified by the agitation rate effects. Since the concentration zones used in this work were calculated at a constant agitation rate of 250 rpm,51 it is expected that there is no exact behavior for experimental concentration when using a time-varying agitation rate trajectory. In Figure 7c (nominal agitation rate trajectory for 1200 μm), the concentration−temperature diagram shows that the system is within the metaestable zone at the beginning of the batch (0−40 min). Then, when the agitation rate is increased, it causes the concentration to reach the unsaturated zone in temperatures around 40 °C, which confirms the results reported by Ni and Liao;24 i.e., the concentration decreases when the agitation rate becomes high. Figure 7d shows the results for the constant agitation rate (case study 3). As shown in this figure, the final D(4,3) calculated by the MOM model (1200 μm) fits appropriately with the experimental data. However, there is a deviation between the MOM model’s prediction and the experimental data at the beginning of the batch due to the fast evaporation of water at the beginning of the process and disturbances into the fluid in the jacket (water steam). After 50 min of batch time, the MOM model fits adequately with the experimental data reaching similar D(4,3). With regards to the FCM, it can be noticed that FCM is not affected by the agitation rate, showing for the experimental case a faster production when compared with the MOM model. Measuring larger D(4,3) and FCM may be associated with crystal sedimentation and inadequate location of the sampling device in the crystallizer. The concentration−temperature diagram for this case shows that, for temperatures below 50 °C, the experimental concentration is located in the unsaturated zone, which allows the dissolution of crystals that have been already formed. Figure 7e shows the results for the robust time-varying agitation rate trajectory of case study 3, where an uncertain region is specified for D(4,3) and FCM. As shown in this figure, the experimental data is located within the uncertain region in D(4,3), which confirms that the proposed methodology for the uncertainty calculation is adequate and allows the behavior of D(4,3) to be predicted in the presence of uncertainty in the system’s kinetic parameters and minimize variability in the batch-to-batch. Remarks. As mentioned above, the agitation rate has no appreciable effect on the FCM production. As mentioned in Section 4.2, the shaded region that covers the FCM variability was calculated using the crystallizer model (eqs 3−15) and the uncertain kinetic parameters obtained from problem 17. Results

from eqs 3−15 only represent an approximation of FCM experimental data as shown in Figure 7b,e. The benefit of using the methodology presented in this work is that the uncertainty in the kinetic parameters has been calculated from experimental data obtained from the pilot-scale process, producing more realistic and practical results than those reported in the literature for batch crystallization.14

5. CONCLUSIONS A study on the agitation rate effects on the production of sugar cane crystals obtained from a pilot-scale process of batch crystallization was presented. The application of the method of moments (MOM) for the PBE solution is adequate since the model predictions are representative of experimental data. The calculation of the kinetic parameters (power-law-type model) is reliable due to use of experimental data from the pilot-scale process under study. The sensitivity analysis performed on both the nucleation rate and the growth rate parameters revealed that the kinetic parameters Kb, p, g, Kg, and q have a significant effect on the D(4,3) and FCM. The methodology developed for the uncertainty analysis while using experimental data was successfully implemented and allowed the identification of a feasible region for the kinetic parameters with significant effect on D(4,3) and FCM. This enabled the specification of robust operating trajectories for the most important variables in batch crystallization (i.e., agitation rate, temperature profile) that can account for product variability in the presence of uncertainty. The agitation rate does not have the same effect over the batch time (the contribution is different at the beginning than at the end of process), since high agitation rate (600−700 rpm) at initial stage of the batch (0−40 min) and a subsequent reduction promotes the crystal growth up to 1000 μm. Furthermore, the opposite agitation rate trajectory must be implemented to reach D(4,3) greater than 1000 μm, i.e., the agitation rate should be low at the beginning of the batch (50− 100 rpm) and increased near the final batch time.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are grateful to the Emerging Leaders in the Americas Program from Canada (ELAP) and the National Council of Science and Technology of Mexico (CONACyTDGEST) for the financial support provided for this research.

■

NOTATION

Capital Latin Letters

A = Total area of heat transfer, cm2 B = Nucleation rate, No. Part./(g·solv·s) Bx = Brix, Bx C = Solute concentration, g·sol/g·sol D(4,3) = Average crystal size in % volume, μm CSD = Crystal size distribution, − ODE = Ordinary differential equation, − PDE = Partial differential equation, − G = Growth rate, cm/min FCM = Formed crystal mass, g MT = Total mass of crystals, g 13192



Article

(12) 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. (13) Mesbah, A.; Landlust, J.; Versteeg, C.; Huesman, A. E. M.; Kramer, H. J. M.; Ludlage, J. H. A; Van den Hof, P. M. J. Model-based optimal control of industrial batch crystallizers. Proceedings of the 20th European Symposium on Computer Aided Process Engineering; Elsevier B.V.: 2010; pp 1563−1568. (14) Nagy, Z. K. Model based robust control approach for batch crystallization product design. Comput. Chem. Eng. 2009, 33, 1685− 1691. (15) Fujiwara, M.; Nagy, Z. K.; Chew, J. W.; Braatz, R. D. Firstprinciples and direct design approaches for the control of pharmaceutical crystallization. J. Process Control 2005, 15, 493−504. (16) Bolaños, E. R. Control and optimization of operating conditions from cooling batch crystallizers. Ph.D. Thesis, I.T. de Celaya, México, 2000. (17) Kim, Y. H.; Lee, K.; Koo, K. K.; Shul, Y. G.; Haam, S. Comparison study of mixing effect on batch cooling crystallization of 3-Nitro-1,2,4-triazol-5-one (NTO) using mechanical stirrer and ultrasound irradiation. Cryst. Res. Technol. 2002, 37, 928−944. (18) Sha, Z.; Louji-Kultanen, M.; Ogawa, K.; Palosaari, S. The effect of mixedness on crystal size distribution in a continuous crystallizer. J. Chem. Eng. Jpn. 1998, 31, 55−60. (19) Kalbasenka, A.; Huesman, A. E. M.; Kramer, H. J. M. Impeller frequency as a process actuator in suspension crystallization of inorganic salts from aqueous solutions. BIWIC 2004: 11th International Workshop on Industrial Crystallization; Hanbat National University: Korea, 2004; pp 135−143. (20) Sander, A.; Kardum, J. P. Pentaerythritol crystallization. Influence of the process conditions on the granulometric properties of crystals. Adv. Powder Technol. 2012, 23, 191−198. (21) Chianese, A.; Kramer, H. J. M. Industrial crystallization process monitoring and control; Wiley-VCH: 2012. (22) Qiu, Y.; Rasmuson, C. Estimation of crystallization kinetics from batch cooling experiments. AIChE J. 1994, 40, 799−812. (23) Saengchan, A.; Kittisupakorn, P.; Paengjuntuek, W.; Arpornwichanop, A. Improvement of batch crystallization control under uncertain kinetic parameters by model predictive control. J. Ind. Eng. Chem. 2011, 17, 430−438. (24) Ni, X.; Liao, A. Effects of mixing, seeding, material of baffles and final temperature on solution crystallization of L-glutamic acid in an oscillatory baffled crystallizer. Chem. Eng. J. 2010, 156, 226−233. (25) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Simultaneous design and control of processes under uncertainty: A robust modelling approach. J. Process Control 2008, 18, 735−752. (26) Ricardez-Sandoval, L. A. Optimal design and control of dynamic systems under uncertainty: A probabilistic approach. Comput. Chem. Eng. 2012, 43, 91−107. (27) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Application of robust control tools to the simultaneous design and control of dynamic systems. Ind. Eng. Chem. Res. 2009, 48, 801−813. (28) Velazquez-Camilo, O.; Bolaños-Reynoso, E.; Lopez-Zamora, L.; Alvarez-Ramirez, J. Experimental evaluation of concentration zone widths in cane sugar crystallization using data and image acquisition. Proceedings of the World Congress on Engineering, London, U.K.; International Association of Engineers, 2010; Vol. 1, pp 709−714. (29) Randolph, A. D.; Larson, M. A. Theory of particulate processes; Academic Press: New York, U.S.A., 1988. (30) Hulburtz, H. M.; Katz, S. Some problems in particle technology: A statistical mechanical formulation. Chem. Eng. Sci. 1964, 19, 555− 574. (31) Gimbun, J.; Nagy, Z. K.; Rielly, D. C. Simultaneous quadrature method of moments for the solution of population balance equations, using a differential algebraic equation framework. Ind. Eng. Chem. Res. 2009, 48, 7798−7812. (32) Dorao, C.; Jakobsen, H. Numerical calculation of the moments of the population balance equation. J. Comput. Appl. Math. 2006, 196, 619−633.

Nr = Agitation rate, rpm P = Pressure, kPa Sr = Supersaturation, − T = Temperature, °C Tj = Jacket temperature, °C U = Global heat transfer coefficient, J/(°C·min·cm2) V = Volume, cm3 Lowercase Latin Letters

b = Exponential of supersaturation in nucleation rate, − g = Exponential of supersaturation in growth rate, − j = Exponential of total mass in nucleation rate, − Kb = Pre-exponential constant in nucleation rate, No. Part./ cm3·min·(g/cm3)j·(rpm)p Kg = Pre-exponential constant in growth rate, cm/min· (rpm)q p = Exponential of agitation rate in nucleation rate, − q = Exponential of agitation rate in growth rate, − n = Population density, No. Part./(μm·g·sol) t = Time, s, min Greek Letters

e = Error, − μ = Distribution moment (volume, surface, length, number), − ρ = Density, g/cm3 Ψs = Uncertain kinetic parameter vector identified from sensitivity analyses, − Ψ̂s = Nominal uncertain kinetic parameter vector identified from sensitivity analyses, −

■

REFERENCES

(1) Jones, A. G. Crystallization process systems; ButterworthHeinemann: 2002. (2) Sarkar, D.; Rohan, S.; Jutan, A. Multi-objetive optimization of seeded batch crystallization process. Chem. Eng. Sci. 2006, 61, 5282− 5295. (3) Hojjati, H.; Sheikhzadeh, M.; Rohani, S. Control of supersaturation in a semibatch antisolvent crystallization process using a fuzzy logic controller. Ind. Eng. Chem. Res. 2007, 46, 1232−1240. (4) Rohani, S.; Horne, S.; Murthy, K. Control of product quality in batch crystallization of pharmaceuticals and fine chemical. Part 1: Design of the crystallization process and the effect of solvent. Org. Process. Res. Dev. 2005, 9, 858−872. (5) Mohd, R.; Baka, A.; Nagy, Z. K.; Saleemi, N. A.; Rielly, D. C. The impact of direct nucleation control on crystal size distribution in pharmaceutical crystallization process. Cryst. Growth Des. 2008, 9, 1378−1384. (6) Mesbah, A.; Nagy, Z. K.; Huesman, M. E.; Kramer, H. J. K.; Van den Hof, J. M. P. Nonlinear model-based control of a semi-industrial batch crystallizer using a population balance modeling framework. IEEE Trans. Control Syst. Technol. 2012, 20, 1188−1201. (7) Lang, Y. D.; Cervantes, A. M.; Biegler, L. T. Dynamic optimization of a batch cooling crystallization process. Ind. Eng. Chem. Res. 1999, 38, 1469−1477. (8) Genk, W. J. Better Growth in batch crystallizers. Chem. Eng. 2000, 90−95. (9) Quintana, P. H.; Bolaños, E. R.; Miranda, B. C.; Salcedo, L. E. Mathematical modeling and kinetic parameter estimation in batch crystallization. AIChE J. 2004, 50, 1407−1417. (10) Bolaños, R. B.; Xaca, O. X.; Á lvarez, J. J. R.; López, L. Z. Effect analysis from dynamic regulation of vacuum pressure in an adiabatic batch crystallizer using data and image acquisition. Ind. Eng. Chem. Res. 2008, 47, 9426−9436. (11) Hu, Q.; Rohani, S.; Jutan, A. Modelling and optimization of seeded batch crystallizers. Comput. Chem. Eng. 2005, 29, 911−918. 13193



Article

(33) Alopaeus, V.; Laakkonen, M.; Aittamaa, J. Solution of population balances with breakage and agglomeration by high-order moment-conserving method of classes. Chem. Eng. Sci. 2006, 61, 6732−6752. (34) Hanks, J. Counting Particles of Cells using IMAQ Vision; Application Note 107; National Instruments, Inc.: Austin, TX, 1997. (35) Mersmann, A. Crystallization Technology Handbook; Marcel Dekker Inc.: New York, U.S.A., 1995. (36) Choong, K. L.; Smith, R. Optimization of batch cooling crystallization. Chem. Eng. Sci. 2004, 59, 313−327. (37) Choong, K. L.; Smith, R. Novel strategies for optimization of batch, semi-batch and heating/cooling evaporative crystallization. Chem. Eng. Sci. 2004, 59, 329−343. (38) Majumder, A.; Nagy, Z. K. Prediction and control of crystal shape distributions in the presence of crystal growth modifiers. Chem. Eng. Sci. 2013, 101, 593−602. (39) Kim, D. Y.; Yang, D. R. A novel method for measurement of crystal growth rate. J. Cryst. Growth 2013, 373, 54−58. (40) Samad, N. A. D. A.; Sin, G.; Gernaey, K. V.; Gani, R. Introducing uncertainty analysis of nucleation and crystal growth models in process analytical technology (PAT) system design of crystallization processes. Eur. J. Pharm. Biopharm. 2013, 85, 911−929. (41) Samad, N. A. D. A.; Sin, G.; Gernaey, K. V.; Gani, R. generic multi-dimensional model-based system for batch cooling crystallization process. Comput. Chem. Eng. 2011, 35, 828−843. (42) Beckman, J. R. Handbook of Chemical Engineering Calculations; McGraw-Hill: New York, 1994. (43) Fernandez, J.; Ballester, L. The heat of crystallization and activity coefficients of sucrose in saturated water solutions. Int. Sugar J. 1971, 73 (866), 40−43. (44) Alvarado, J. S. Optimization of vapor system operation in the batch crystallization process, using a response surface method. MSc Thesis. I.T. de Orizaba, México, 2009. (45) Yakimento, O. A. Engineering computations and modeling in MATLAB/Simulink; AIAA Educations Series; American Institute of Aeronautics and Astronautics: 2011. (46) Ricardez-Sandoval, L. A. Current challenges in the design and control of multiscale systems. Can. J. Chem. Eng. 2011, 89, 1324−1341. (47) Evans, R. D.; Ricardez-Sandoval, L. A. Multi-scenario modelling of uncertainty in stochastic chemical systems. J. Comput. Phys. 2014, 273, 374−392. (48) Rasoulian, S.; Ricardez-Sandoval, L. A. Uncertainty analysis and robust optimization of multiscale process systems with application to epitaxial thin film growth. Chem. Eng. Sci. 2014, 116, 590−600. (49) Luna, P. A. E. Adaptive control strategy for a vacuum batch crystallizer. MSc Thesis, I.T. de Orizaba, México, 2014. (50) Akrap, M.; Kuzmanic, N.; Kardum, P. J. Effect of mixing on the crystal size distribution of borax decahydrate in batch cooling crystallizer. J. Cryst. Growth 2010, 312, 3603−3608. (51) Antonio, A. A. Experimental determination of the concentration zones of sugar cane molasses for experimental implementation through virtual instrumentation. MSc Thesis, I.T. de Orizaba, México, 2011.

13194


Dynamic Modeling and Optimization of Batch Crystallization of Sugar

Recommend Documents