Mathematical modelling and analysis of an industrial scale evaporative crystallizer producing lactose monohydrate

Journal of Food Engineering 154 (2015) 49–57

Contents lists available at ScienceDirect

Journal of Food Engineering journal homepage: www.elsevier.com/locate/jfoodeng

Mathematical modelling and analysis of an industrial scale evaporative crystallizer producing lactose monohydrate Shailesh Agrawal a,⇑, Anthony (Tony) Paterson a, Jeremy McLeod b, James Jones a, John Bronlund a a b

School of Engineering and Advanced Technology, Massey University, Palmerston North, New Zealand Hilmar Cheese Company, Hilmar, CA, USA

a r t i c l e

i n f o

Article history: Received 4 July 2014 Received in revised form 22 November 2014 Accepted 22 December 2014 Available online 10 January 2015 Keywords: Lactose mono-hydrate Crystallization Industrial Model Simulation Parameter estimation Crystallizer dynamics MSPMR Size dependent growth

a b s t r a c t Lactose is industrially produced by a series of concentration steps to supersaturate the solution allowing crystallization to take place. This papers deals with the operation and analysis of a forced circulation evaporative crystallizer producing pharmaceutical lactose. A mathematical model of the operation of the crystallizer was developed. The model comprises of the mass balance, population balance, nucleation kinetics, growth kinetics and lactose solubility and mutarotation rate expressions. The steady state operation of the crystallizer was analysed from the classical mixed suspension mixed product removal (MSMPR) model perspective. It was found that the crystallizer behaviour deviated from ideal MSMPR behaviour with a steep upward curvature in the population density curve for crystal sizes below 40 lm. A size dependent growth (SDG) rate model was used to model the deviation from ideal behaviour. The parameters of the SDG model were estimated by fitting it to the measured population density curve from the crystallizer. The crystallization kinetic parameters like growth and nucleation rate constants were then estimated by fitting the simulated trends to the measured date. The model showed that the average crystal size follows a transient decaying dynamics before settling down to a constant value. This was reflected in the actual data collected from the crystallizer. The transient dynamics disappeared from the model results when a simple size independent growth model was used. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Lactose is generally produced by concentration of whey permeate followed by a series of batch cooling crystallisers. Mathematical models have been developed for cooling batch lactose crystallization by numerous researchers (McLeod, 2007; Mimouni et al., 2009; Vu et al., 2006; Wong et al., 2010). During the permeate concentration step, the supersaturation must remain below that where spontaneous nucleation occurs. This limits the concentration that can be achieved and thus the yield (Paterson, 2009).Evaporative crystallisation does not have this limitation, as the continuous crystallization that occurs, exhausts the supersaturation created due to the water removal by evaporation. This allows the yield to be increased. However, evaporative crystallization is a process of continuous nucleation and growth that leads to a broader particle size distribution than that in the batch crystallization with a single nucleation event. Understanding the dynamics and crystallisation processes occurring in the evaporative ⇑ Corresponding author at: Department of Chemical Engineering, Visvesvaraya National Institute of Technology, Nagpur 440010, India. Tel.: +91 9545378898. E-mail address: [email protected] (S. Agrawal). http://dx.doi.org/10.1016/j.jfoodeng.2014.12.025 0260-8774/Ó 2015 Elsevier Ltd. All rights reserved.

crystallisers will enable better control of the crystal size distribution (CSD) coming out of the evaporative crystalliser. This work presents, to the best of authors’ knowledge, the first detailed study of an industrial scale forced circulation evaporative crystallizer for lactose monohydrate. 2. Mathematical model development In the development of this model the following assumptions have been made: (i) The crystallizer operates under ideal mixed suspension mixed product removal (MSMPR) conditions. It is well mixed and the temperature, crystal suspension and supersaturation are uniform throughout the crystallizer. The product removal is isokinetic and no classification of the crystals occurs. (ii) Lactose crystals are spherical. (iii) No agglomeration or crystal breakage occurs. This assumption was supported by the microscopic images of the lactose crystals which showed intact, single crystals coming out of the crystallizer.

50

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

Nomenclature Ca Cb Cas CLS C Cc Bs Bp D[4,3] F G Ge Km L M Q T b g j kevap kn,p kn,s kg m m0 m1

dissolved a-lactose concentration [kg (kg of solution)1] dissolved b-lactose concentration [kg (kg of solution)1] equilibrium solubility of a-lactose [kg (kg of water)1] equilibrium solubility of lactose [kg (kg of water)1] total dissolved lactose concentration [kg (kg of solution)1] crystal content of slurry [kg (kg of slurry)1] secondary nucleation rate [# min1 (kg of slurry)1] primary nucleation rate [# min1 (kg of slurry)1] volumetric weighted mean diameter (lm) correction factor for a-LMH solubility growth rate (lm min1)) effective size independent growth rate in the SDG model of Mydlarz (lm min1) mutarotation rate constant bin size/crystal size (lm) mass of the slurry present in the crystallizer (kg) mass flow rate (kg min1) crystallizer temperature (°C) secondary nucleation rate order growth rate order exponential dependence of secondary nucleation on crystal content evaporation rate (kg min1) primary nucleation rate constant [# min1 (kg of slurry)1((kg of a-LMH) (kg of water)1)b] empirical secondary nucleation rate constant [# min1 (kg of crystal)1((kg of a-LMH) (kg of water)1)b] growth rate constant [lm min1 ((kg of a-LMH)(kg of water)1)g] mass of single crystal (kg) 0th momentof CSD [# (kg of slurry)1] 1st moment of CSD [m (kg of slurry)1]

(iv) The orders for growth and nucleation kinetics from the literature can be applied to the industrial crystallizer under study i.e. they are scale independent. (v) The fluctuations and disturbances to flow rates and temperatures common during industrial operation are neglected.

2.1. Mass balance A schematic representation of an evaporative crystallizer is shown in Fig. 1. Q(kg min1), u[(kg of solution) (kg of slurry)1], C[(kg of dissolved lactose) (kg of crystal free slurry)1] represents the mass flow rate, voidage and dissolved lactose concentration respectively. Subscripts f and o denotes the feed and the discharge streams to and from the crystallizer. M(kg) is the mass of the slurry present in the crystallizer at any given time. kevap(kg min1) is the evaporation rate. The overall lactose solute balance in the crystallizer is given by:

m2 m3 m4 n(L) np t

2nd moment of CSD [m2 (kg of slurry)1] 3rd moment of CSD [m3 (kg of slurry)1] 4th moment of CSD [m4 (kg of slurry)1] population density of crystal size L [# (kg of slurry)1 (lm)1] primary nucleation rate order time (min1)

Greek symbols q density (kg m3) u voidage [(kg of solution) (kg of slurry)1] 1u crystal content of slurry [(kg of crystalline lactose) (kg of slurry)1] s residence time (min) Subscripts f inlet/feed stream o outlet/product stream c crystal s lactose solution S slurry w water CSD crystal size distribution LMH lactose monohydrate MOL method of lines MSMPR mixed suspension mixed product removal crystallizer ODEs ordinary differential equations PDEs partial differential equations GRD growth rate dispersion SDG size dependent growth TS total solids

The industrial crystallizer being modelled operates in two distinct modes. During the start-up, lactose solution is fed into the crystallizer vessel until a targeted volume is achieved. The steam flow to the calandria (heat exchanger) is then initiated, providing energy to concentrate the solution. The feed flow is controlled to maintain a constant volume in the crystallizer. Fresh lactose solution is continuously fed to the crystallizer to make up for the reduction in volume due to evaporation. The start-up phase can therefore be considered as a semi-batch or fed batch operation. This semi-batch concentration process continues till the desired slurry density is achieved. Then the discharge valve is opened and slurry is pumped forward to the next stage in the process. At this point the feed flow is set to maintain a constant level in the crystallizer, allowing for evaporation and the discharge flow. This marks the start of the continuous operation. Eq. (1) for a semi-batch system under no discharge (Qo = 0) and crystal free feed (uf = 1) conditions, becomes:

d d ½M uC þ ½0:95Mð1  uÞ ¼ Q f uf C f  Q o uo C o þ 0:95Q f ð1  uf Þ dt dt  0:95Q o ð1  uo Þ ð1Þ where the terms on the left hand side represent the rate of accumulation of the dissolved and crystalline lactose in the crystallizer. The first two terms on the right hand side give the feed and discharge rate of dissolved lactose and the last two terms gives the feed and discharge rate of crystalline lactose. The crystalline lactose terms are multiplied by 0.95 to account for the water of crystallization.

kevap Qf,φf,Cf

M,φ,C

Qo,φo,Co

Fig. 1. A generalized schematic of an evaporative crystallizer.

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

d d ½M uC þ ½0:95Mð1  uÞ ¼ Q f C f dt dt

ð2Þ

Eq. (2) on differentiating and rearranging takes the form:

dC M½0:95  C ¼ dt

du dt

 ½uC þ 0:95ð1  uÞ dM þ Q f Cf dt Mu

ð3Þ

Solving Eq. (3) gives the dissolved lactose concentration profile for the semi-batch process. The rate of change of crystallizer mass during this process is given by Eq. (4), where kevap is the mass evaporation rate.

dM ¼ Q i  kev ap dt

ð4Þ

The feed rate of solution was not directly measured in the industrial crystallizer under study. However, since the crystallizer was operated at constant volume and during batch operation no discharge was occurring, the volumetric feed rate was equal to the volumetric rate of water removal, as shown in Eq. (5), where qw is the density of water and qs is the density of feed solution.

Qi

qs

¼

kev ap

qw

ð5Þ

The evaporator operates as a single effect with no vapour recompression allowing the evaporation rate to be taken as the steam mass flow rate fed to the calandria neglecting heat losses. In the continuous mode, M is a constant as the crystallizer operates at constant volume and density. For the continuous mode, under the conditions of crystal free feed (uf = 1) and well mixed assumption (Co = C, uo = u), Eq. (1) on rearrangement becomes:

dC Q f C f C 0:95ð1  uÞ ðC  0:95Þ du ¼    dt M u s su u dt

ð6Þ

where s is the mean residence time (min) given by ratio of mass in the crystallizer to the discharge flow rate, expressed as:

s¼

M Qo

ð7Þ

This leaves u and ddtu, the slurry voidage terms to be defined. The voidage is defined by the crystal content of the slurry and is therefore governed by the number and size of the crystals present in the slurry. This leads to the population balance equation (PBE). 2.2. Population balance The unsteady state PBE for a continuous constant volume, isothermal, well mixed crystallizer with no agglomeration and breakage is given by Randolph and Larson (1988) as

@nðL; tÞ @nðLÞ nðLÞ ¼ G  : @t @L s

ð8Þ

The term on the left hand side (LHS) of Eq. (8) gives the rate of change of crystal number density where n(L) is the population density of crystal size L[# (kg of slurry)1(lm)1]. L and t refer to the size (lm) and time (min1) domain respectively. The first term on the right hand side (RHS) gives the net change in crystal number density due to growth in and out of the size interval dL, where G (lm min1) is the size independent growth rate given by dL/dt. Size dependent growth rate is discussed later while analyzing the industrial crystallizer data. The second term on the RHS of Eq. (8) gives the number of crystals within the size interval dL that enter and leave the system with the inflow and the outflow of the slurry. The symbol s is the residence time of the crystallizer (min). The above population balance equation can be solved subjected to the following boundary and initial conditions:

51

Bp;s G nð1; tÞ ¼ 0

ð10Þ

nðL; 0Þ ¼ ni ðLÞ

ð11Þ

nðL0 ; tÞ ¼

ð9Þ

Eq. (9) represents the boundary condition at the lower end of the size distribution (L0). The production of the nuclei at the smallest size (L0) occurs (at all time t) by nucleation given by Bp,s/G, where B is the nucleation rate or the birth rate of nuclei [# min1(kg of slurry)1]and G is the growth rate (lm min1). Subscripts p and s represents primary and secondary nucleation mechanisms, respectively. No direct nucleation occurs in the higher size ranges. The boundary condition for the upper end of the size distribution, Eq. (10), is chosen such that the number of crystals at that size is zero. It is imperative for the mass balance to hold that the crystal size does not outgrow this maximum limit. Eq. (11) gives the size distribution of the crystals present in the crystallizer, at time t = 0. Method of lines (MOL) is a commonly used technique (Abbas and Romagnoli, 2007; Griffin et al., 2010; Quintana-Hernández et al., 2004) for solving the population balance PDEs. In MOL, the PDE is discretized along the linear size domain, resulting in the formation of a set of ODEs which can be solved numerically. The discretization can be done using various finite difference schemes. The following set of ODEs has been derived using a second order finite difference scheme in length as they show less oscillations than higher order approximation while providing satisfactory numerical efficiency and accuracy (Griffin et al., 2010).

9 dnðL0 Þ 3nðL0 Þ þ 4nðL1 Þ  nðL2 Þ > > ¼ > > dL 2 DL > > = dnðLi Þ nðLiþ1 Þ  nðLi1 Þ ¼ ; i ¼ 1; 2; 3 . . . p  1 > dL 2 DL > > > dnðLp Þ 4nðLp Þ  7nðLp1 Þ þ 4nðLp2 Þ  nðLp3 Þ > > ; ¼ dL 2 DL

ð12Þ

where DL = Li–Li1 and p is the number of discretized intervals and p + 1 is the number of grid points. The first grid point n(L0) corresponds to the boundary condition in Eq. (3). The PBE scheme discussed above requires growth and nucleation kinetic models which are discussed in Section 2.3. 2.3. Lactose crystallization kinetics During the start-up process of an evaporative crystallizer, in the absence of any seeding, crystal formation occurs by primary nucleation. After the start up, the system contains a slurry of dissolved and crystalline lactose and operates at a lower supersaturation. During this stage secondary nucleation becomes the dominant mechanism of nucleation. The primary nucleation, secondary nucleation and growth kinetics of alpha lactose monohydrate (a-LMH) can be expressed as

np Ca  C as ; 1C  b Ca  C as ; Bs ¼ kn;s C jc 1C

Bp ¼ kn;p



ð13Þ ð14Þ

and

G ¼ kg



Ca  C as 1C

g ð15Þ

As secondary nucleation occurs in the presence of solute crystals, the concentration of crystals is included in the rate expression with the order of dependence, j. np, b and g represent the order of dependence of primary nucleation, secondary nucleation and growth kinetics, respectively on supersaturation (the bracketed term) and were fixed based on the review of existing literature.

52

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

There are numerous studies investigating growth rate of lactose monohydrate. For this work, growth rate order (g) of 1.3 as determined by McLeod (2007) was used. A limited literature exists on the secondary nucleation kinetics of lactose. Analysis by Liang et al. (1991) was the most robust. From their study, the value of secondary nucleation rate order (b) was fixed as 1.5. In absence of any other study, data from Mcleod (2007) was used to determine the value of primary nucleation rate order (np) as 2.5. kn,p [m min1((kg of a-LMH) (kg of water)1)1.3], kn,s[# min1 (kg of slurry)1((kg of a-LMH) (kg of water)1)2.5] and kg[# min1 (kg of slurry)1((kg of a-LMH) (kg of water)1)1.5], are the rate constants and these were estimated by fitting the model predicted values to the measured industrial crystallizer data. Agrawal (2012) has a detailed review of the lactose crystallization kinetics. The parameters required for the calculation of supersaturation driving force in alpha-lactose monohydrate crystallization are discussed next. The equilibrium a-lactose solubility, Cas [kg (kg of water)1] in the presence of b-lactose is given by Visser (1982) as

C as

 C  C LS  FK m 1C  C LS ¼ ; 1 þ Km

ð16Þ

where F is a temperature dependent correction factor suggested by Visser (1982) to account for the suppression of solubility of a-lactose in presence of b-lactose. Its dependence on temperature, T (°C), can be given by the following exponential relationship

F ¼ 0:0187exp0:02367T

ð17Þ

1

CLS [kg (kg of water) ] is the lactose equilibrium solubility and can be estimated at a given temperature by the following expression (McLeod, 2007)

C LS ¼

10:9109exp0:02804T : 100

ð18Þ

The total dissolved lactose concentration (C) is the sum of a and b lactose concentrations,

C ¼ Ca þ Cb:

ð19Þ

The equilibrium mutarotation rate constant, Km is defined as

Km ¼

Cb : Ca

ð20Þ

Combining Eqs. (19) and (20), the a-LMH concentration can be calculated from the total lactose concentration by

Ca ¼

C : ð1 þ K m Þ

ð21Þ

McLeod (2007) fitted the data by Roetman and Buma (1974) to give the following equation to predict Km as a function of temperature.

K m ¼ 0:0024T þ 1:6353

ð22Þ

2.4. Slurry properties: crystal content, mean crystal size and density The voidage of the slurry can be calculated as

u ¼ 1  qc m3 :

ð23Þ

The second term on the RHS of Eq. (23) gives the crystal content of the slurry. m3 is the third moment of the number density distribution and is a measure of the volume of the total number of crystals present in the slurry [(m3 of crystal) (m3of slurry)1]. It is defined as

m3 ¼

Z

1

L3 ndL:

ð24Þ

0

Differentiating Eq. (23) gives

du dm3 ¼ qc ; dt dt where dm3/dt is given by Randolph and Larson (1988)

ð25Þ

dm3 m3 ¼ 3Gm2  : dt s

ð26Þ

m2 in Eq. (26) is the second moment of the density distribution and defined as

m2 ¼

Z

1

L2 ndL

ð27Þ

0

The volume weighted mean diameter D[4,3](Randolph and Larson(1988)) is given by

D½4; 3 ¼

m4 : m3

ð28Þ

m4 in the above equation is the fourth moment of the density distribution,

m4 ¼

Z

1

L4 ndL

ð29Þ

0

The crystal content Cc [units] of the slurry is calculated as,

C c ¼ m3 qc :

ð30Þ

Total solids of the slurry can be calculated as,

TS ¼ uC þ 0:95C c :

ð31Þ

A bench top study found that the density of the slurry remained unchanged irrespective of the state of solute (dissolved or crystalline). Hence, a relationship between the dissolved solids and the solution density can be used to relate total solids to the slurry density. Due to a lack of data on density of lactose solution at higher concentrations, density data of sucrose solution from Bubnik et al. (1995) at the crystallizer temperature was used. The density of sucrose solution matched very closely (a maximum error of 1%) with the available density data for lactose solutions (Buma (1980)).

qs ¼ 983:2e0:0042ðTSÞ

ð32Þ

The slurry density in itself does not present a complete picture regarding the crystal yield as it is a function of total solids, which in turn depends on both dissolved as well as crystalline solids. Hence, the crystal concentration needs to be evaluated separately in order to estimate the production rate of crystals. 3. MSMPR analysis of industrial crystallizer It is possible to ascertain crystallization kinetics in an evaporative crystallizer using the steady state mixed suspension mixed product removal (MSMPR) criterion. Under these conditions a plot of the natural log of the population density and crystal size results in a straight line. From the slope and the intercept of the line, growth and nucleation rate respectively can then be calculated. For steady state operation to be assumed, Randolph and Larson (1988) state that approximately eight residence times are required. The residence time in the crystallizer was recorded directly from the control room HMI (Human Machine Interface) consoles (Fig. 2A). To ensure steady state for the crystallizer being studied, sampling was started at least 15 h after the crystallizer was switched to continuous mode. Once this state was achieved in the crystallizer samples of the slurry were collected every 30– 40 min for 210 min. Where a stoppage occurred sampling was abandoned because the steady state was disturbed. Three non-stop sampling events under steady state were achieved over a period of a month. Measurements were also conducted during the transient regime of the crystallizer immediately after it was moved from the semi-batch to continuous mode for approximately 8 h. At each sampling time slurry was collected in sealed plastic 125 ml sampling containers. The sample was split and a portion transferred into a separate sample container, cooled under running

53

Residence me (hours)

1.5

Trial 1

Trial 2

Dissolved solids (kg per kg of crystal free slurry)

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

Trial 3

1.3 1.1 0.9 0.7 0.5 0

200

0.43

Trial 1, MA

Trial 2, MA

0.42

Trial 3, MA

Trial1, HRM

Trial 2, HRM

Trial 3, HRM

0.41 0.4 0

Trial 2

Trial 3

0.72 0.7 0.68 100

150

200

Trial 1

0.45

0.4 50

100

150

Time (minutes)

Time (minutes)

(C)

(D)

30 mins

6

80 mins

Trial 1

Trial 2

200

250

Trial 3

380 360

D[4,3] (μm)

7

Trial 3

0.5

400

135 mins

4 165 mins

3

Trial 2

0.55

0

0 min

5

200

0.6

250

9 8

100

(B)

0.74

Volume %

0.44

(A)

0.76

50

0.45

Time (minutes)

Trial 1

0

0.46

Time (minutes)

Crystalline solids (kg of crystals per kg of slurry)

Total solids (kg per kg of slurry)

0.78

100

0.47

340 320 300 280 260

205 mins

2

240

1

220 200

0 1

10

100

1000

0

100

Size (μm)

Time (minutes)

(E)

(F)

200

Fig. 2. Profiles of the key parameters during steady state operation of an industrial crystallizer (A) residence time, (B) dissolved solids (MA: moisture analyzer and HRM: handheld refractometer), (C) total solids, (D) crystalline solids (calculated from mass balance using total solids and dissolved solids), (E) volume weighted mean diameter and (F) volumetric crystal size distribution for trial 1.

water and stored for total solids determination. The cooling was done to reduce evaporative losses. The remaining slurry was allowed to stand for two minutes to allow the crystals and liquid to separate. The liquid layer was siphoned off using a 10 ml syringe and 1–1.5 g of liquor was filtered through a 0.44 lm filter on to the pan of the halogen moisture analyzer (HB43, Mettler Toledo). The total dissolved solids were also measured using a handheld digital refractometer (Atago, PAL-1). The moisture analyzer readings were considered more reliable because of the dependence of refractive index on temperature and the composition of the dissolved solids. A comparison of the two values obtained by the two techniques is shown in Fig. 2B. Following this, the total solids of the sample was determined (Fig. 2C). The slurry was well mixed with a help of a spatula and 0.8–1.2 g of sample was smeared over the moisture

metre pan. The analyzer quickly heats upto 155 °C and it takes 6–10 min for the analysis to be complete, depending on the mass of sample. Knowing the total solids and dissolved solids, the crystalline content of the slurry was determined by mass balance (Fig. 2D). The slurry sample was also used for crystal size determination using a Malvern Mastersizer (2000 S). The slurry was diluted with slightly supersaturated lactose solution at 25 °C (20% w/w) to reduce the viscosity of the slurry so that a homogenous suspension of the crystal mass could be obtained. The Malvern cell was flushed twice with water and once with the dispersant lactose solution before each analysis. Fresh dispersant (lactose solution) was then introduced and stirring started and the background measurement was allowed to stabilize. The slurry sample was shaken to ensure

54

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

uniform suspension of the crystals and introduced into the Malvern cell using a transfer pipette. Since the particle concentration of the slurry was high, 5–10 drops of sample were sufficient to achieve the desired obscuration levels in the Malvern cell (10– 20%). The measured CSD for trial 1 is shown in Fig. 2E. The volume weighted mean diameter for all the three measurement trials are plotted in Fig. 2F. The measured parameters can be considered fairly stable over the three hour period. Using these results it can be concluded that the crystallizer generally operates with the following steady state parameters in continuous mode: residence time: 1 h; dissolved solids (mass fraction): 0.45–0.46; total solids (mass fraction): 0.72–0.73; crystalline content (mass fraction): 0.49–0.52; D[4,3]: 300–340 lm. Having established the steady state operation, in order to estimate the crystallization kinetics, the volumetric distribution needs to be converted to a number density distribution. This can be done using the following formula (modified from Tait et al., 2009))

nðLÞ ¼

6C c v

ð33Þ

100pqc L3 ðDLÞ1018

where v, is the mass percent occupied by crystal of size L (Malvern output); L the geometric mean size of the two consecutive size bins used in Malvern output (lm); DL the difference between the two consecutive size bins used in Malvern output (lm); qc the density of a-lactose monohydrate crystal (kg m3) and Cc the crystal content in terms of weight fraction of the slurry [(kg of crystal) (kgof slurry)1]. The plot of ln(n(L)) versus L for the CSD measured under steady state for trial 2 is shown in Fig. 3. The plot is not a straight line. It has a strong upward curvature at about 40 lm, making it difficult to calculate accurate crystallization kinetics from the MSMPR model. A linear extrapolation leads to an under estimation of nucleation rate while over estimating the growth rates (Liang et al., 1991). This deviation from the MSMPR model can be explained by growth rate dispersion (GRD) or size dependent growth rate (SDG) (Liang et al., 1991). GRD dispersion refers to a distribution of growth rates experienced by a crystal population; due to either different conditions they experience in a crystallizer or due to some inherent structural conditions that differ from crystal to crystal (Randolph and Larson, 1988). SDG refers to the dependence of growth rate on the crystal size with growth increasing with size. SDG and GRD are not linked to a specific mechanism; rather they are just a phenomenological classification of non-ideal growth behaviour which can be caused

25

20

ln (n)

15

10

5

0 0

200

400

600

800

1000

1200

Size (μm) Fig. 3. Steady state population density curve obtained from the industrial crystallizer at steady state and the fitted growth model.

by a number of mechanisms such as difference in dislocation density when a spiral growth mechanism is dominant, the presence of strain within the crystal lattice, Ostwald ripening, size dependent diffusion constant and attrition of crystals (Menon,2006). A comprehensive review by Garside (1980) discusses SDG and GRD. Both GRD and SDG can lead to the same crystal size distribution and thus, it is not possible to distinguish between the two phenomena by MSMPR size analysis alone (Randolph and Larson, 1988). There are numerous GRD and SDG models in the literature. Mydlarz (1996) fitted many of these various GRD and SDG models to the population density curve, concluding that the SDG models do a better job of predicting the population density curves and suggesting further development of GRD models as there is experimental evidence to support GRD phenomenon. In this study, a SDG model was used to describe the upward curvature in the population density curve because of the ease of applying this model. 3.1. Size dependent growth rate model The size dependent model used was proposed by Mydlarz (1993) as follows:

GðLÞ ¼ Ge f1  exp½aðL þ cÞg

ð34Þ

When substituted in the population balance equation for an ideal MSMPR system, the population density equation takes the following form:

nðLÞ ¼ no expðaLÞ

 1ðasGe Þ=ðasGe Þ exp½aðL þ cÞ expðacÞ  1

ð35Þ

where Ge is the effective size independent growth rate (lm min1); no the ‘zero-size’ population density [# lm1 (kg of slurry)1] and a and c are the model parameters. DataFit (Oakdale Engineering, trial version 9) was used to fit Eq. (35) to the plant data to estimate the parameters for the SDG model. The fitted values of the parameters a, c and Ge for the Mydlarz model was 0.0059, 0.0091 lm and 1.82 lm min1, respectively. The fitted model is shown in Fig. 3. 4. Simulation and parameter estimation The equations for the crystallizer were solved in MATLAB (R2007) using the ode45 routine. Size dependent growth model, Eq. (34), was used to account for the deviation from MSMPR behaviour. The effective size independent growth rate (Ge) in the SDG model was expressed similar to the normal growth rate expression in Eq. (15). The semi-batch and continuous mode equations were solved separately with the end conditions of the semi-batch passed as the initial conditions for the continuous mode. The semi-batch and the continuous mode were simulated for 60 and 3940 min of operation, respectively. No seeding was used in the studied crystallization process and hence initial condition, Eq. (11), becomes n(L,0) = 0. For the boundary conditions it was assumed that all nuclei were formed at a size of 1 micron. This assumption is supported by data presented by Shi et al. (1989). The population balance was solved using the method of lines (MOL). In MOL, the size domain is discretized into p number of grid points. The accuracy of the solution of the resulting ODEs increases with the increase in the number of grid points. However, a greater number of grid points also increases the computational time. Hence a balance needs to be achieved between the desired accuracy and the computational time. Simulations were conducted at 100, 200, 500, 1000, 2000 and 4000 grid points. The maximum crystal size was chosen as 1000 lm based on the actual crystal size distribution measured from the plant. Hence the grid sizes

55

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57 Table 1 Kinetic parameters estimated from size dependent and size independent growth models. Kinetic parameter ? Model ;

kg [m min1((kg of a-LMH) (kg of water)1)1.3]

kn,p [# min1 (kg of slurry)1((kg of a-LMH) (kg of water)1)2.5]

kn,s [# min1 (kg of slurry)1 ((kg of a-LMH) (kg of water)1)1.5]

SDG SIG Literature value

42 32 7

8E8 3E8 3.1E8

4E10 1.3E8 1E10

Fig. 4. Fitted and measured profiles for the key parameter of an industrial forced circulation evaporative lactose crystallizer. (A) Dissolved lactose, (B) total solids, (C) crystal concentration (D) volume weighted mean diameter and (E) steady state population density distribution (at time t = 4000 min). The dotted line through the measured data points in (D) is drawn as a visual guide to compare with the model prediction. The broken and continuous lines represent the model prediction for size dependent and size independent growth models, respectively. The symbols represent the measured values. The fitted parameters are growth, primary nucleation and secondary nucleation rate constants.

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

corresponded to 10, 5, 2, 1, 0.5 and 0.25 lm, respectively for the number of grid points used. It was observed that below 1 lm there was no significant improvement in the simulation results. A grid size of 1 lm and Lmax = 1000 lm was used for all following simulations. The understanding gained by the sensitivity analysis (Agrawal, 2012) was used to optimize the various kinetic parameters to match the model predictions with the measured steady state data. The supersaturation exponents for growth, primary and secondary nucleation processes were kept constant at 1.3, 1.5 and 2.5 respectively and were assumed to be independent of scale of operation. The rate constants were estimated by trial and error to achieve the best fit by the following considerations. The total solids value depends on the evaporation rate and hence evaporation rate was the first parameter to be set based on total solids measurements. The secondary nucleation rate constant drove the particle size and the dissolved solids concentration (supersaturation). It was thus optimized to approximately predict the measured steady state mean particle size (325 lm) and dissolved lactose (0.435 mass fraction). The growth rate constant was then tuned to get a closer fit to the measured points. The first measurements during the transient regime provided the crystal size at the end of the batch process (100–130 lm) and were used to fix the primary nucleation rate constant. The fitted parameters are shown in Table 1. The simulations showed crystal size exhibiting transient behaviour before reaching a steady state. To validate this sampling was conducted immediately after the crystallizer was switched to continuous mode for 8 h. Total solids, dissolved lactose and CSD were determined for the collected samples. The measured data points are shown along with the simulation results in Fig. 4. The discontinuous line denotes the profile during the batch process. No sampling during the batch stage was possible and hence this could not be validated. Sampling for 8 h was insufficient to cover the entire transient stage as the model shows that the crystal size requires approximately 25 h to stabilize. The other parameters are stable long before this. The decaying transient dynamics of the crystal size cannot be validated by direct measurements as it was not possible to sample continuously for 25 h. However, an alternative way to confirm whether such dynamics are exhibited by the crystallizer was established. Normal plant operation procedures require the operators to sample the dried lactose crystals and perform a sieve test using an 80 and 200 mesh sieve every couple of hours. The sieve test data was collected after three different start-ups. The results are shown in Fig. 5. The increase in the mass fraction of crystals retained on the 80 mesh (170 lm) sieve indicates an increase in the mean crystal size and vice a versa. It is clear that the crystal size shows oscillatory dynamics. The first cycle of increasing particle size, reaching a maxima and then dropping to a minimum size is easily distinguishable. Although exact crystal sizes were not measured, the trends support the model findings. Such dynamics have also been reported by numerous previous simulation and experimental studies of evaporative industrial crystallizers (Bermingham et al., 1998; Kramer et al., 1996, 2000). These dynamics also probably explain the reason for the smaller mean particle size measurement in trial 1 as they may have been measured during the time when the crystal size was passing through the lower point of the oscillation cycle. The measured crystal size cycle is out of sync with the model predictions but follows the trend accurately. The final steady state crystal density distribution was well predicted by the model using the optimized parameters as is shown in Fig. 4E. It is to be noted that the parameters were not optimized to fit the population density curve. A good fit proves the validity of the model.

90 80

% retained at 80 mesh screen

56

70 60 50 40 30 20 10 0 0

500

1000

1500

Time (minutes) Fig. 5. Dry crystal sieve analysis representative of the CSD dynamics of the evaporative crystallizer.

Such dynamics can be caused by numerous factors such as fines destruction, classified product removal, size dependent growth rate, kinetics of nucleation rate and external modulation of parameters (Lakatos et al., 2007). Thus it is difficult, if not impossible, to exactly figure out the source of instability (Randolph and Larson, 1988). In the present system, no fines destruction or classified product removal occurred, therefore size dependent growth rate or growth rate dispersion could have caused the dynamic behaviour. Simulations with size independent growth (SIG) rate were conducted; the results are shown (by broken lines) in Fig. 4. The kinetic parameters were adjusted to make the simulation results match as closely as possible with the measured data. The fitted growth and nucleation rate constant values along with the values reported in the corresponding literature from which the orders for growth and nucleation were fixed in the model are shown in Table 1. Considering the stochastic nature of crystallization and the magnitude of difference in the scales at which these were determined, the fitted values are in reasonable agreement with the literature values It can be seen that the decaying transient disappears when the growth becomes size independent and the model failed to predict the fluctuations in the crystal size observed in the plant. Therefore, SDG captures the plant dynamics more precisely than SIG. An expected straight line is obtained on a population density versus crystal size plot for SIG as shown in Fig. 4E. It is to be noted that selection of the SDG model does not suggest that the a-LMH crystals exhibit size dependent growth rate. More research is warranted on this aspect. The SDG model was selected because of the better fit obtained and its simpler mathematical form compared to the GRD models (Agrawal, 2012). 5. Conclusion A mathematical model for an industrial evaporative crystallizer was formulated and simulated using MATLAB. The population density curve for the industrial crystallizer showed strong upward curvature for crystal sizes below 40 lm. Size dependent growth (SDG) rate was incorporated into the model to account for this deviation. Lactose growth and nucleation rate constants were estimated by fitting the model to the measured plant parameters. It was found that the crystal mean size from the crystallizer showed dampening oscillatory dynamics after the start-up and then gradually settled down to a steady state. The developed mathematical model was able to capture this behaviour successfully. Size independent

S. Agrawal et al. / Journal of Food Engineering 154 (2015) 49–57

growth model did not give rise to any such dynamics showing that size dependent growth or growth rate dispersion needs to be used when modelling these systems. Acknowledgement The authors will like to thank Foundation for Research Science and Technology, the Government of New Zealand and Fonterra Cooperative Group for funding the project. References Abbas, A., Romagnoli, J.A., 2007. Multiscale modeling, simulation and validation of batch cooling crystallization. Sep. Purif. Technol. 53 (2), 153–163. Agrawal, S.G., 2012. Evaporative crystallization of alpha-lactose monohydrate. (PhD Thesis), Massey University, Palmerston North. Bermingham, S.K., Kramer, H.J.M., van Rosmalen, G.M., 1998. Towards on-scale crystalliser design using compartmental models. Comput. Chem. Eng. 22 (Supplement 1), S355–S362. Bubnik, Z., Kadlec, P., Urban, D., Bruhns, M., 1995. Sugar Technologists Manual. Chemical and Physical Data for Sugar Manufacturers and Users. Bartens Pub. Co, Berlin. Buma, T.J., 1980. Viscosity and density of concentrated lactose solutions and of concentrated cheese whey. Netherland Milk Dairy J 34, 65–68. Garside, J., Davey, R.J., 1980. Invited review secondary contact nucleation: kinetics, growth and scale-up. Chem. Eng. Commun. 4 (4), 393–424. Griffin, D.W., Mellichamp, D.A., Doherty, M.F., 2010. Reducing the mean size of API crystals by continuous manufacturing with product classification and recycle. Chem. Eng. Sci. 65 (21), 5770–5780. Kramer, H.J.M., Dijkstra, J.W., Neumann, A.M., Ó Meadhra, R., van Rosmalen, G.M., 1996. Modelling of industrial crystallizers, a compartmental approach using a dynamic flow-sheeting tool. J. Cryst. Growth 166 (1–4), 1084–1088. Kramer, H.J.M., Dijkstra, J.W., Verheijen, P.J.T., Van Rosmalen, G.M., 2000. Modeling of industrial crystallizers for control and design purposes. Powder Technol. 108 (2–3), 185–191. Lakatos, B.l.G., Sapundzhiev, T.J., Garside, J., 2007. Stability and dynamics of isothermal CMSMPR crystallizers. Chem. Eng. Sci. 62 (16), 4348–4364.

57

Liang, B., Shi, Y., Hartel, R.W., 1991. Growth rate dispersion effects on lactose crystal size distributions from a continuous cooling crystallizer. J. Food Sci. 56 (3), 848– 854. http://dx.doi.org/10.1111/j.1365-2621.1991.tb05397.x. McLeod, J., 2007. Nucleation and growth of alpha lactose monohydrate. (PhD Thesis), Massey University, Palmerston North, New Zealand. Menon, A.R., 2006. A task based design procedure and modelling approaches for industrial crystallization processes. (PhD Thesis), Delft University of Technology, Netherland. Mimouni, A., Schuck, P., Bouhallab, S., 2009. Isothermal batch crystallization of alpha-lactose: a kinetic model combining mutarotation, nucleation and growth steps. Int. Dairy J. 19 (3), 129–136. http://dx.doi.org/10.1016/ j.idairyj.2008.09.006. Mydlarz, J., 1996. Modelling of growth rate for MSMPR crystallizer data. Cryst. Res. Technol. 31 (5), 541–565. http://dx.doi.org/10.1002/crat.2170310502. Mydlarz, J., Jones, A.G., 1993. On the estimation of size-dependent crystal growth rate functions in MSMPR crystallizers. Chem. Eng. J. Biochem. Eng. J. 53 (2), 125–135. Paterson, A.H.J., 2009. Production and uses of lactose. In: Fox, P.F., McSweeney, P.L.H. (Eds.), Advanced Dairy Chemistry, vol. 3. Springer Science + Business Media LLC, New York, pp. 105–121. Quintana-Hernández, P., Bolaños-Reynoso, E., Miranda-Castro, B., Salcedo-Estrada, L., 2004. Mathematical modeling and kinetic parameter estimation in batch crystallization. AIChE J. 50 (7), 1407–1417. Randolph, A.D., Larson, M.A., 1988. Theory of Particulate Processes. Academic Press, San Diego. Roetman, K., Buma, T.J., 1974. Temperature dependence of the equilibrium b/a ratio of lactose in aqueous solution. Netherland Milk Dairy J 28, 155–165. Shi, Y., Hartel, R.W., Liang, B., 1989. Formation and growth phenomena of lactose nuclei under contact nucleation conditions. J. Dairy Sci. 72 (11), 2906–2915. Tait, S., White, E.T., Litster, J.D., 2009. A study on nucleation for protein crystallization in mixed vessels. Cryst. Growth Des. 9 (5), 2198–2206. Visser, R.A., 1982. Supersaturation of alpha-lactose in aqueous solutions in mutarotation equilibrium. Netherlands Milk Dairy J 36, 89–101. Vu, T.T.L., Durham, R.J., Hourigan, J.A., Sleigh, R.W., 2006. Dynamic modelling optimisation and control of lactose crystallisations: comparison of process alternatives. Sep. Purif. Technol. 48 (2), 159–166. http://dx.doi.org/10.1016/ j.seppur.2005.07.015. Wong, S.Y., Bund, R.K., Connelly, R.K., Hartel, R.W., 2010. Modeling the crystallization kinetic rates of lactose via Artificial Neural Network. Cryst. Growth Des. 10 (6), 2620–2628. http://dx.doi.org/10.1021/cg100122y.