Comparison of FEM and DPB Numerical ... - ACS Publications

ARTICLE pubs.acs.org/IECR

Comparison of FEM and DPB Numerical Methodologies for Dynamic Modeling of Isothermal Batch Gibbsite Crystallization Andrey V. Bekker, Tian S. Li, and Iztok Livk* CSIRO Process Science and Engineering and Light Metals Flagship, Parker CRC for Integrated Hydrometallurgy Solutions, Waterford WA 6152, Australia ABSTRACT: A population balance equation based dynamic gibbsite crystallization model, incorporating crystal growth, nucleation, and agglomeration, was solved using two different numerical techniques, namely, a fully implicit finite element method (FEM) and an explicit discretized population balance (DPB) numerical scheme. Unlike previous crystallization modeling approaches, the agglomeration model in the FEM framework was formulated based on the Safronov agglomeration equation and its partial differential equation (PDE) approximation [Bekker, A. V.; Livk, I. Agglomeration process modeling based on a PDE approximation of the Safronov agglomeration equation. Ind. Eng. Chem. Res. 2011, in press. The FEM numerical solution is implemented using a fully implicit Newton’s method Galerkin finite element algorithm, which applies automatic Gear-type timestep and nonuniform adaptive mesh strategies to help optimal solution convergence. The dynamic FEM and DPB model predictions of isothermal gibbsite crystallization, using estimated kinetic parameters, are validated against experimental crystallization data that were generated under process conditions relevant to Bayer gibbsite crystallization. Additional modeling results obtained for a modified set of kinetic parameters, prolonged crystallization times, and more complex crystal size distributions show that the novel FEM-based crystallization modeling framework offers a more accurate and computationally efficient model solution than that based on the DPB approach.

1. INTRODUCTION Processes that involve solid particles, liquid droplets, or gas bubbles suspended in a continuous phase represent a class of particulate systems that are very common across a range of processing and manufacturing industries. They include various unit operations ranging from gravity separations to crystallization and precipitation. Complex behavior of such systems requires relatively sophisticated tools to enable the development of accurate process models needed for effective process design, scale-up, optimization, and model-based control. In this work, we focus on modeling crystal size distributions (CSDs) in gibbsite crystallization, an operation that represents a crucial step in the production of alumina and aluminum by the Bayer process. Commonly, in gibbsite (aluminum trihydroxide) crystallization a number of different mechanisms affect the evolution of the CSD, including crystal growth, nucleation, and agglomeration. The population balance equation (PBE), introduced for the description of crystallization processes by Hulburt and Katz2 and Randolph and Larson,3 provides an elegant mathematical framework for coupling of crystal size distribution and crystallization kinetics. As PBE-based models of realistic processes are rarely amenable to analytical solutions, they have to be solved numerically in most cases. This has become considerably more attainable in recent years due to the advent of more effective numerical solution algorithms and ever-increasing computing power. There have been four main approaches to the solution of PBE models: (i) finite difference methods, (ii) finite element methods (FEM), (iii) pivotal schemes,4 and (iv) the method of moments. A discretized population balance (DPB) represents one of the finite difference solution methods. The most widely used form of Published 2011 by the American Chemical Society

the DPB, which also incorporates the description of crystal agglomeration, was developed by Hounslow et al.5 Using a different approach, Gelbard and Seinfeld6 developed a collocation method—one of the first finite element method based solution algorithms. Their approach has been subsequently extended by a number of other authors.7,8 Roussos and coauthors9 used a fixed order polynomial basis and showed that stability was improved when compared to the orthogonal collocation and discrete population balance solutions. An alternative solution of the PBE, as described by Randolph and Larson,3 is based on the moments of particle number distribution and can considerably reduce the complexity of the model solution. A further simplification is achieved by implementing the quadrature method of moments introduced by McGraw.10 In the case where the complete CSD is required, it has to be reconstructed from the moments through additional calculations. The DPB and FEM methods, however, resolve a complete CSD directly. The focus of this work is on applying a novel FEM framework for modeling gibbsite crystallization1,11 and comparing the new solution to a well-established DPB model solution.5,12 The novel FEM numerical approach consists of two modules: (i) nucleation and growth module11 and (ii) agglomeration module.1 The two modules are consolidated in this work to provide a unified framework for modeling gibbsite crystallization, with all three crystallization mechanisms taking place simultaneously. Unlike Received: September 20, 2010 Accepted: December 23, 2010 Revised: November 29, 2010 Published: February 28, 2011 3994

dx.doi.org/10.1021/ie101942s | Ind. Eng. Chem. Res. 2011, 50, 3994–4002

Industrial & Engineering Chemistry Research

ARTICLE

previous crystallization models, the current FEM-based model has been formulated using the Safronov agglomeration equation13 and its new partial differential equation (PDE) approximation, as derived for the gibbsite crystallization by Bekker and Livk.1,14 They showed that the use of a new agglomeration model led to a fully implicit FEM framework, which combined Newton’s method and an implicit-in-time Galerkin finite element algorithm with fixed order polynomial basis, nonuniform adaptive mesh, and automatic Gear-type adaptive time step. As there is no analytical solution available for a dynamic population balance crystallization model with simultaneous nucleation, growth, and agglomeration, the novel FEM model solution is compared to a DPB model solution. DPB and FEM model predictions of the gibbsite crystallization process, including evolutions of the CSD and aluminate concentration, are compared in this work to assess the accuracy and efficiency of the novel FEM method for solving crystallization models. Model predictions are validated against experimental crystallization data to verify that the new model can be applied for describing the behavior of real-world crystallization systems.

Furthermore, as discussed by Bekker and Livk,14 dynamic batch crystallization modeling can be applied to approximate the behavior of a series of continuous crystallizers, as employed in the production of alumina using the Bayer process. A dynamic 1-D PBE model of isothermal well-mixed batch gibbsite crystallization that couples the number density of crystal size distribution and the supersaturation balance is given as ∂nðt; vÞ ∂ þ ðGðt; vÞ nðt; vÞÞ ¼ Bu ðtÞδðv-vmin Þ ∂t ∂v ∂nðt; vÞ þ ; v∈½0;vmax ; t∈½0;tend ∂t a Z v ∂nðt; vÞ 1 ¼ β0 ðtÞ nðt; v - uÞ nðt; uÞ du ∂t 2 0 a Z vmax nðt; uÞ du -β0 ðtÞ nðt; vÞ 0

Z ð1-

2. PB MODELING OF GIBBSITE CRYSTALLIZATION The particulate phase can be characterized by a number of different properties such as, size, shape, surface area, mass, and chemical composition. In crystallization systems, commonly, crystal size and the total number of crystals change with time or location. As the system involves a large number of crystals of different sizes, a continuous PBE formulation is used to describe the evolution of crystal size distribution. In a general form,3 the population balance equation can be stated in the particulate phase space as ∂n þ r 3 wj n ¼ B - D ∂t

unðt; uÞ duÞkM 0

dσðtÞ dt

¼ - vmin Bu ðtÞ- 3kV 1 = 3 GL ðtÞ μ2 = 3 ðtÞ; Z μ2 = 3 ðtÞ ¼

vmax

u2 = 3 nðt; uÞ du

ð2Þ

0

Gðt; vÞ ¼

ð1Þ

where n is the number density distribution, t is time, and wj is the velocity in the direction of the jth internal coordinate such as crystal size. The terms B and D represent the birth and death functions, which can account for the crystallization processes, such as nucleation and agglomeration. For the purpose of this work the internal coordinate in eq 1 is taken as crystal diameter or volume. In gibbsite crystallization the crystal size distribution of the crystalline product is typically affected by a number of crystallization mechanisms. The crystal growth mechanism enlarges particle diameter via the deposition of gibbsite from caustic aluminate solution on existing crystals. At the same time, nucleation increases the total number of crystals in the system through the generation of new entities. Agglomeration, on the other hand, reduces the total number of crystals by the cementation of a number of small crystals into a larger single crystal. All these three crystallization mechanisms affect the total surface area available, with only nucleation and agglomeration affecting the total number of crystals. Crystal breakage would also increase the number of crystals in the system, but this mechanism is not considered in this work, as typically breakage of crystals does not play an important role in gibbsite crystallization.15 2.1. A 1-D Dynamic PBE Model of Batch Gibbsite Crystallization. Although a continuous crystallization process is commonly used in the production of bulk chemicals such as gibbsite, a batch laboratory-scale crystallizer can be used to perform studies of industrial crystallization systems on a smaller scale.

vmax

dv ¼ 3kV 1 = 3 v2 = 3 GL ðtÞ; GL ðtÞ ¼ Kg σg ðtÞ; dt Bu ðtÞ ¼ Kn σn ðtÞ; β0 ðtÞ ¼ Ka σa ðtÞ

nðt; vÞjv ¼ 0 ¼ 0; nðt; vÞjt ¼ 0 ¼

∂nðt; vÞ j ¼0 ∂v vmax

0; v < vmin ; n0 ðvÞ; v g vmin

σðtÞjt ¼ 0 ¼ σ0

where n is the number density, t is time, v and u are the crystal volume, subscript “a” denotes the agglomeration term, and σ is relative supersaturation. Note that in the above model the supersaturation balance equation is strictly valid only for an isothermal case. n0 is the initial number density of seed. G is the volume-based crystal growth rate, with GL = dL/dt the linear crystal growth rate that is assumed independent of crystal size.16 μ2/3 is the moment of a volume-based CSD that is proportional to the total surface area of crystals. Bu is the secondary nucleation rate, and kM is the conversion factor defined as 2M[Al(OH)3]Aeq/(FsM[Al2O3]). The description of other model parameters and their values are given in Tables 1 and 2. β0 is the size-independent agglomeration kernel15 required for the description of the agglomeration rate. Next, the general population balance model of eq 2 is adjusted separately for the DPB and the FEM solution framework. 2.2. A DPB Gibbsite Crystallization Model. A discretized population balance approach5 has been used widely to solve complex population balance based models involving nucleation, agglomeration, and crystal growth. The method approximates a continuous size coordinate, as used in the original PDE (eq 1), with a finite number of discretized crystal size intervals, resulting 3995

dx.doi.org/10.1021/ie101942s |Ind. Eng. Chem. Res. 2011, 50, 3994–4002


ARTICLE

Table 1. Parameters Used in the Model of Isothermal Gibbsite Crystallization parameter

meaning

Lmin vmin = kVLmin

3

value

dimension

minimal crystal size

5.024 10-6

m

minimal crystal volume

6.64 10-17

m3

-6

Lmax

maximal crystal size

200 10

vmax = kVLmax3

maximal crystal volume

4.19 10-12

m3 s

m

tend

final process time

6000

kV

volume shape factor

π/6

kM

conversion factor

4.46 10-2

vn V0

nuclei dispersion analytical seed exponential factor

vmin/10 1.19 10-14

m3 m3

N0

seed total crystal number

6.85 1011

no. m-3

M0

seed total crystal mass

20

kg m-3

σ0

initial relative supersaturation

0.83

Fs

gibbsite density

2420

kg m-3

M[Al(OH)3]

molecular weight of aluminum trihydroxide

77.99

g mol-1

M[Al2O3]

molecular weight of aluminum oxide

101.94

g mol-1

Aeq

equilibrium solubility of Al2O3 at experimental conditions

74.15

kg m-3

Table 2. Estimated Values of Gibbsite Crystallization Kinetic Parameters parameter

meaning

value

dimension

Kg

growth rate coefficient

8.2 10-9

m s-1

Kn

nucleation rate coefficient

1.24 108

m-2s-1

K0

agglomeration kernel

5.28 10-15

m3 s-1

g n

growth rate order nucleation rate order

3.3 5

-

a

agglomeration rate order

6

-

technique has been used before to model gibbsite crystallization, as reported by Li et al.,12,19,20 demonstrating its ability to accurately model the behavior of this crystallization system. The DPB approach is therefore chosen as a reference model to which a novel FEM-based gibbsite crystallization model will be compared. 2.3. FEM Gibbsite Crystallization Model. A FEM-based gibbsite crystallization model will be built here by consolidating two separate modules of a continuous PBE model of gibbsite crystallization that have been developed recently by Bekker and Livk.1,11,14,21 Using particle volume as an internal coordinate, the FEM-based model is stated as

in a finite number of ordinary differential equations. A DPB crystallization model, given as dNi dNi dNi dNi ¼ þ þ ð3Þ dt dt g dt n dt a was implemented to resolve the evolution of crystal size distribution in a batch gibbsite crystallization system by incorporating the contributions of crystal growth, nucleation, and agglomeration to the rate of change of the number of crystals in the ith size interval, as indicated in eq 3. The discretized growth term, (dNi/dt)g, where N denotes the number of crystals per volume of slurry, was evaluated using a two-term discretized growth model as described by Hounslow et al.5 The discretized nucleation source term, (dNi/dt)n, was implemented as ( Bu , i ¼ 1 dNi ¼ ð4Þ 0, otherwise dt n Finally, the discretized agglomeration term, (dNi/dt)a, was evaluated using a discretized agglomeration formulation proposed by Hounslow et al.5 and later extended by Litster et al.17 and Wynn18 to allow for a finer discretization of the crystal size coordinate. A geometric size discretization with the volume ratio of viþ1/vi = 21/q and q = 6 was chosen for the DPB crystallization simulations run in this work. The resulting system of ordinary differential equations (ODEs) was solved using a fourth order Runge-Kutta adaptive time step integration method. The DPB

∂nðt; vÞ ∂ þ 3kV 1 = 3 GL ðv2 = 3 nðt; vÞÞ ∂t ∂v 2 Bu v-vmin ∂nðt; vÞ ¼ pffiffiffi exp ð5Þ þ ∂t vn π vn a Z v ∂nðt; vÞ ∂ ¼ -β0 ½nðv; tÞ unðu; tÞ du ∂t ∂v 0 a Z vmax -β0 nðv; tÞ nðu; tÞ du v

where Bu is the secondary nucleation rate that is multiplied by a narrow Gaussian-shape source rate function with an average nuclei size of vmin and nuclei size dispersion of vn. This function is used to approximate the original delta function used traditionally to constrain the size range of the nucleation rate.22 GL is the sizeindependent linear crystal growth rate, and β0 is the sizeindependent agglomeration kernel, with the rate expressions as those presented in eq 2. The last term on the right-hand side in the second line of eq 5 represents the contribution of crystal agglomeration. The formulation of the agglomeration rate is based on the Safronov agglomeration equation, which was adapted for use within the FEM-based modeling framework by Bekker and Livk.1 To make it fit into the implicit FEM framework, the original Safronov agglomeration equation was transformed into a PDE agglomeration 3996



ARTICLE

form as ∂nðt; vÞ ∂ ¼ -β0 ðnðt; vÞ μ1 ðt; vÞÞ - β0 nðt; vÞðμ0 ðt;vmax Þ ∂t ∂v a -μ0 ðt; vÞÞ ∂μ0 ðt; vÞ ¼ nðt; vÞ; μk ðt; vÞ ¼ ∂v

Z

v

uk nðt; uÞ du

ð6Þ

0

∂μ1 ðt; vÞ ¼ vnðt; vÞ ∂v In the above PDE agglomeration model, n is the number density. μ0, μ1, and μk are zero, first, and kth cumulative moments of the number density function, defined within a given size range. The mass balance equation and the equation for the total surface area of crystals, eq 2, were transformed to a form amenable to the implicit numerical scheme as ð1-μ1 ðt;vmax ÞÞkM

dσðtÞ ¼ -vmin Bu ðtÞ - 3kV 1 = 3 GL ðtÞ μ2 = 3 ðt;vmax Þ dt

ð7Þ ∂μ2 = 3 ðt; vÞ ¼ v2 = 3 nðt; vÞ ∂v The resulting 1-D dynamic isothermal batch crystallization model, eqs 5-7, was solved using a general adaptive fully implicit Newton’s method based Galerkin finite element algorithm. The applied algorithm used automatic Gear-type time-step and nonuniform adaptive mesh strategies, which helped achieve optimal solution convergence. This numerical strategy leads to a fully implicit (i.e., time advanced) numerical algorithm with a band-diagonal sparse Jacobian matrix. The implicit time discretization of the mass balance equation (eq 7) also provides a sparse band-diagonal Jacobian matrix. The details of the implementation of the FEM-based numerical methodology to solve the gibbsite crystallization model are given by Bekker and Livk.1,11 The FEM solution of the gibbsite crystallization model is compared to the DPB solution obtained using exactly the same crystallization kinetics in both models. To validate the ability of the new FEM modeling framework to describe a real-world crystallization system, model predictions are compared to experimental gibbsite crystallization data, which are described in section 3.

3. GIBBSITE CRYSTALLIZATION EXPERIMENTAL DATA 3.1. Isothermal Batch Gibbsite Crystallization Experiment. Gibbsite is the polymorph of Al(OH)3 that is crystallized from supersaturated caustic aluminate solutions in the commercial Bayer process used for extracting alumina from bauxite. Experimental conditions, such as temperature and solution concentration, applied in a laboratory isothermal gibbsite crystallization experiment were typical of the agglomeration section of gibbsite crystallization circuits as implemented in alumina refineries. Laboratory crystallization experiments were conducted in a 4 L jacketed stainless-steel crystallizer fitted with a draft tube and an 86 mm axial flow Lightnin A310 impeller. A schematic of the crystallizer and other details of the experimental setup can be

found in Li et al.20 The isothermal batch crystallization experiment reported here was carried out at 80 °C, with an agitation rate of 480 rpm and a volume mean shear rate of 459 s-1. Synthetic caustic aluminate solution used in the experiment was prepared by digesting technical grade Al(OH)3 in an aqueous sodium hydroxide solution, with sodium carbonate of 40 g L-1 added to the aluminate solution. The aluminate solution was pressure filtrated and allowed to cool and equilibrate at temperature, T, of 80 °C. This resulted in an aluminate concentration, A, of 140 g L-1 Al2O3 and a caustic concentration, C, of 200 g L-1 Na2CO3. The total alkali concentration was 240 g L-1 Na2CO3. The equilibrium solubility23 of the synthetic aluminate solution at those conditions was Aeq = 74.15 g L-1 Al2O3, which gives a relative supersaturation, σ = A/Aeq - 1, of 0.89. In the model, a slightly lower initial relative supersaturation value of 0.83, as measured at the end of the induction time, is used. Gibbsite was being crystallized from caustic aluminate solution by seeding with 15 g L-1 Alcoa C31 gibbsite crystals initially. The laboratory crystallizer is well-mixed24 in both liquor and solid phase. Samples were taken periodically for the determination of the suspension density, liquor composition, and CSD. The CSD at each sampling time was determined using a Coulter Counter Multisizer, which measures the number-based crystal size distribution. 3.2. Estimation of Gibbsite Crystallization Kinetics. A differential method12,25 was used in this work for estimating gibbsite crystallization kinetics from batch experimental data. The estimation is based on the CSDs and solution supersaturation values measured at different sampling times during a 100 min crystallization experiment. The initial induction period was not taken into account. For each time instant, the zero and third moments, which are equivalent to the total number of crystals and total crystal volume, respectively, were estimated from measured CSD data. They were used together with the number of crystals in the first size interval to estimate the three crystallization kinetic rates for nucleation, Bu, crystal growth, GL, and agglomeration, β0. The estimated kinetic rates were correlated with the corresponding values of supersaturation in order to formulate a simple power-law relationship between the three crystallization kinetics and solution supersaturation. The values estimated for the kinetic parameters are presented in Table 2. A comprehensive study on the estimation of gibbsite crystallization kinetics was reported previously by Li et al.26

4. MODEL PREDICTIONS A 1-D dynamic isothermal well-mixed batch gibbsite crystallization model of eq 2 was solved with the estimated crystallization kinetics using the DPB and FEM methodologies, as stated by eqs 3 and 4 and eqs 5-7, respectively. In each case, the experimental initial CSD and initial supersaturation of 0.83 defined the required starting values for the model. The initial CSD resulted in a total number of crystals, Ntotal, of 6.85 1011 no. m-3 and the total crystal mass, Mtotal, of 20 kg m-3. Using estimated gibbsite crystallization kinetics, four different cases were solved using the DPB and FEM numerical approaches, including (i) all three crystallization mechanisms active, (ii) nucleation mechanism only, (iii) agglomeration mechanism only, and (iv) crystal growth mechanism only. While the first case enables a direct comparison of the model predictions to experimental data, the other three cases provide additional insight about the consistency of the results obtained by the 3997



ARTICLE

Figure 1. Experimental data and model predictions for (a) CSDs at final time and (b) supersaturation evolutions.

Figure 2. DPB and FEM predicted CSDs for (a) nucleation only, (b) agglomeration only, and (c) crystal growth only cases at constant supersaturation and 100 min residence time.

two different modeling approaches. The results from the first case are shown in Figure 1, with the results from the other three cases shown in Figure 2. To convert the FEM volume number density results to the number of crystals in each size interval, used in the DPB model and shown in Figure 1, the following relationships were used: Z viþ1 nðuÞ du Ni ¼ vi

nL ¼ 3kV 1 = 3 v2 = 3 nðvÞ

ð8Þ

Ni is the number of crystals in the ith size interval with a volume of vi - viþ1, where viþ1 = vi 21/6, vi = kV Li3, and L1 = Lmin = 5.024 μm. nL is the number density function in the L crystal size space. 4.1. Comparison between DPB/FEM Model Predictions and Experimental Data. DPB- and FEM-based numerical predictions of the gibbsite batch experiment are presented in terms of the final CSD and the evolution of solution concentration in Figure 1. Considering the simplicity of the kinetic rate equations implemented in modeling, the CSDs predicted by the two models agree reasonably well with the experimental CSD at the final time, as shown in Figure 1a. The FEM predicted CSD can better match the narrow peak of the experimental CSD than the DPB model prediction, which exhibits a lower maximum due 3998



ARTICLE

Table 3. Differences between the Final and Initial Values of Total Number and Total Mass of Crystals As Predicted by the Two Models and Measured in the Experiment DPB Figure

Bu

GL

β0

1

Knσn

Kgσg

Kaσa

2a

Knσ0n

0

0

2b 2c

0 0

0 Kgσ0g

Kaσ0a 0

-3

ΔNtotal, no. m

-3.07 1011

FEM -3

ΔMtotal, kg m

-3

ΔNtotal, no. m

experiment -3

ΔMtotal, kg m

ΔNtotal, no. m-3

ΔMtotal, kg m-3

-2.91 1011

33.7

31.2

-3.04 1011

30.8

11

0.0

11

2.95 10

0.0

-

-

-5.35 1011 0.32 1011

0.0 113.3

-5.35 1011 0.01 1011

0.0 102.9

-

-

2.95 10

Figure 3. DPB and FEM model predictions of (a) CSD and (b) supersaturation evolution, for a 20-h crystallization run.

to a higher degree of numerical diffusion. Overall, the shape of the FEM predicted CSD, shown in Figure 1a, matches the experimental data more closely than the CSD predicted by the DPB model. At the same time, both methods predict a very similar supersaturation evolution, as shown in Figure 1b. FEM and DPB predictions of the total number and the total mass of crystals are also very similar, exhibiting up to 5% deviation from the experimental data. Values of the deviations in the measured and model predicted total number and total mass of crystals are collated in Table 3 for all four cases mentioned above. 4.2. Numerical Tests for Cases with Single-Mechanism Kinetics. To gain additional insight into the difference between the DPB and FEM numerical solutions, three other hypothetical cases, with only one crystallization mechanism taking place at a time, were solved using the two models at constant supersaturation conditions. In the nucleation only case, the crystal growth rate and agglomeration kernel parameters were set to zero; therefore, Kg = Ka = 0 in eq 2. The nucleation mechanism only affects the first size interval, which is reflected in an increased number of crystals in that interval at the final time, as presented in Figure 2a. For the nucleation only case, the number of crystals in the second interval is expected to remain unchanged. Although the two modeling approaches give very similar results for the first size interval, it is evident that the FEM nucleation only solution leads to a slight overprediction of the numbers in the second interval when compared to the result predicted by the DPB model. This is a consequence of the Gaussian approximation (in eq 5) of the nucleation rate delta function (in eq 2) as implemented in the FEM model. However, the difference between the total number of crystals predicted by the two

models is less than 0.1%, as shown in Table 3, for the nucleation only case. Both models accurately predict a zero increase in the total mass of crystals. For the agglomeration only case, with Kn = Kg = 0, similar CSDs were predicted by the two modeling techniques, as shown in Figure 2b. A small difference between the FEM and DPB predicted CSDs is a result of the Safronov agglomeration equation, implemented within the FEM methodology, which is only an approximation of the agglomeration process modeled by the Smoluchowski agglomeration equation and incorporated within the DPB model. Both models agree well in the prediction of the total number of crystals and total mass of crystals, as shown in Table 3. In the growth only case, with Kn = Ka = 0, the FEM and DPB solutions, presented in Figure 2c, were obtained for a constant crystal growth rate. When the size-independent crystal growth is the only active mechanism, the CSD translates horizontally to the right without changing its shape. For this case, as shown in Figure 2c, the DPB number density function shows a higher level of numerical diffusion than that predicted by the FEM model. This results in a lower maximum of the CSD in the DPB case. The FEM model, which utilizes a nonuniform adaptive mesh technique, maintains a more consistent shape of the initial number density function, exactly as expected for a crystal growth only case. Compared to the FEM solution, numerical diffusion in the DPB case leads to around 5% overestimation of the total number of crystals and 8% overestimation of the total mass of crystals, as indicated in the last row of Table 3. The results presented in this section clearly show that the FEM solution is associated with a lesser degree of numerical diffusion in the crystal growth only case than that of the DPB model, resulting in a better preservation of the CSD shape and its 3999



ARTICLE

maxima, which is reflected in more accurate prediction of the total number of crystals and total mass of crystals. 4.3. Prolonged Crystallization Residence Time. To further demonstrate the difference between the FEM- and DPB-based solutions, the same gibbsite crystallization process was simulated for a residence time of 20 h. Results from this simulation, as presented in Figure 3, show increased difference in the CSD predictions by the DPB and FEM models. The difference between the two models increased due to the accumulation of the numerical diffusion error in the DPB case. In contrast, the evolutions of relative supersaturation predicted by the two models are similar. The agreement between the final total crystal mass and the total number of crystals predicted by the two models is better than 5%, as shown in Table 4. The FEM model, however, required only one-third the CPU time of the DPB model. 4.4. Increased Nucleation Rate and Reduced Minimum Crystal Size. Additional numerical tests were conducted to compare the DPB and FEM numerical algorithms for the cases with much higher nucleation rates and smaller minimum crystal sizes. Crystal growth and agglomeration kinetic parameters were kept the same as those estimated from the experimental data in Table 2. A gibbsite nucleation rate 10 times higher than that estimated from the gibbsite experimental data was used to simulate gibbsite crystallization with a residence time of 100 min. In these cases, the experimental initial CSD was approximated by a negative exponential function: n0(v) = (N0/V0) exp(-v/V0). The total seed mass and total number of crystals of the approximated initial CSD are the same as those of the experimental CSD. Numerical experiments were run for two different minimum sizes, Lmin. Model predicted final CSDs are shown in Figure 4. At the initial time, the high nucleation rate creates an instant CSD maximum at the minimum crystal size, which then slowly propagates to a larger crystal size range due to growth and agglomeration. The CSD maximum develops primarily due to a Table 4. Total Number and Total Mass of Crystals, and CPU Time Required To Obtain Model Solution on a 2.4 GHz Processor, for a Residence Time of 20 h model

Ntotal, m-3

Mtotal, kg m-3

CPU time, s

DPB

3.78 10

11

92.8

90

FEM

3.63 1011

92.0

30

difference between the initial value of the number density and the ratio of Bu/GL at the minimum crystal size. This phenomenon is discussed in more detail by Bekker and Livk.11 At the same time, the main peak of the initial CSD shifts to the right and reduces in magnitude due to crystal growth and agglomeration. Presented in Figure 4b is a case where the original CSD is extended to a minimum crystal size of 0.5 μm. Comparing the magnitudes of the nucleation peaks predicted for two different Lmin values, as presented in Figure 4, it is evident that the two peak heights are the same and therefore independent of the nuclei size. Reducing the nuclei size from Lmin = 5 μm to Lmin = 0.5 μm simply shifts the nucleation CSD peak position to the left by 4.5 μm, while the position of the right-hand-side CSD peak remains the same. If the agglomeration mechanism is eliminated by setting β0 = 0, a growth and nucleation case study is obtained, which can be validated against an analytical solution. For a given supersaturation profile, this case results in a number density, which can be expressed as a function of time and size11 as nL ðL, tÞ 8 Rt Rt g > > ΔL ¼ 0 GL ðtÞ dt ¼ Kg 0 σ ðtÞ dt > < B ðt -tÞ K u end n ¼ σ n - g ðtend -tÞ, at L ¼ Lmin þΔL if L e Lmin þ ΔL ¼ G ðt -tÞ K L end g > > > : nL ðL þ ΔL, tÞ, if L > Lmin þ ΔL

ð9Þ Analytical and numerical CSDs, evaluated for the case of growth and nucleation only, are presented in Figure 5a. Simultaneously with the CSDs, the supersaturation profiles were calculated by the DPB and FEM numerical algorithms, as shown in Figure 5b. The FEM supersaturation profile was taken into account in the evaluation of the analytical CSD, which was estimated using eq 9. Both model predicted CSDs match the analytical solution relatively closely. However, the FEM model solution tracks the analytical solution much more closely than the DPB model. A superior performance of the FEM model is particularly notable at the points where the number density function exhibits large gradients. The DPB model has difficulty in matching the sharp maximum of the number density function, predicting a much lower amplitude than the analytical solution. The lower amplitude in the DPB model case is mainly due to numerical diffusion associated with the fixed DPB grid. Unlike the DPB numerical

Figure 4. DPB and FEM predicted CSDs for Bu = 10Knσn, for a minimum crystal size of (a) 5 μm and (b) 0.5 μm after 100 min of crystallization. 4000



ARTICLE

Figure 5. DPB and FEM model predictions for (a) final CSD and (b) evolution of supersaturation, for the case with Bu = 10Knσn, β0 = 0, and the measured crystal growth rate, for a minimum size of 0.5 μm.

algorithm, the FEM numerical approach applies an adaptive nonuniform size discretization scheme, which can optimally refine the mesh at regions with high gradients in order to achieve the required accuracy. Modeling results presented in this section show that DPB and FEM approaches lead to similar model predictions, but the FEM predictions matched the experimental data of a 100 min residence time experiment more closely than the DPB predictions. For the extended parameter range and for longer residence times, the FEM algorithm was less diffusive, could more closely track the complex crystal size distribution, and was computationally more efficient than the DPB-based solution algorithm.

5. CONCLUSIONS A dynamic process model of isothermal batch gibbsite crystallization, formulated using a PBE mathematical framework, was solved for Bayer process conditions implementing two different numerical approaches, namely, a fully implicit Newton’s method based Galerkin FEM algorithm and an explicit DPB numerical approach. Using the values for gibbsite crystallization kinetic rates based on laboratory experimental data, numerical results from both DPB and FEM models showed good agreement with the measured CSD and supersaturation data. All three main crystallization mechanisms, nucleation, crystal growth, and agglomeration, were shown to affect the shape of the CSD of the gibbsite crystallization product significantly. A more detailed comparison between the FEM and DPB numerical solutions suggests the following conclusions: (i) The CSD predicted by the FEM numerical solution for the case with complete gibbsite crystallization kinetics showed less numerical diffusion and matched the experimental data more closely than the CSD predicted by the DPB model. In terms of the total number and the total mass of crystals, for this case the FEM and DPB numerical results agreed within a 5% deviation for up to 20 h of crystallization time. At the same time, the FEM numerical algorithm was found to be around 3 times more efficient than the DPB-based algorithm. (ii) Simulation results, obtained for a hypothetical case with the crystal growth kinetics only, clearly showed that, unlike the DPB solution, the FEM model solution was able to preserve a consistent shape of the CSD, as expected for this case based on analytical reasoning.

(iii) The nuclei crystal size assumed by the model does not affect the amplitude of the CSDs at its left boundary. The maximum of the nucleation CSD peak position is shifted by the difference in the nuclei size. (iv) For more complex multipeak CSDs, the FEM model solution was found to match the analytical solution significantly more closely than the DPB model solution. The superior performance of the FEM based solution in that case is mainly the result of a nonuniform adaptive size discretization strategy used by the FEM numerical algorithm.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ REFERENCES (1) Bekker, A. V.; Livk, I. Agglomeration process modeling based on a PDE approximation of the Safronov agglomeration equation. Ind. Eng. Chem. Res. 2011, DOI: 10.1021/ie101933r. (2) Hulburt, H. M.; Katz, S. Some problems in particle technology. A statistical mechanical formulation. Chem. Eng. Sci. 1964, 19, 555–574. (3) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization, 2nd ed.; Academic Press: Toronto, 1988. (4) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization—I. A Fixed Pivot Technique. Chem. Eng. Sci. 1996, 51, 1311–1332. (5) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. A discretized population balance for nucleation, growth and aggregation. AIChE J. 1988, 34 (11), 1821–1832. (6) Gelbard, F.; Seinfeld, J. H. Numerical solution of the dynamic equation for particulate systems. J. Comput. Phys. 1978, 28, 357–375. (7) Nicmanis, M.; Hounslow, M. J. Finite-element method for steady-state population balance equations. AIChE J. 1998, 44, 2258– 2272. (8) Wulkow, M.; Gerstlauer, A.; Nieken, U. Modeling and simulation of crystallization processes using parsival. Chem. Eng. Sci. 2001, 56, 2575–2588. (9) Roussos, A. I.; Alexopoulos, A. H.; Kiparissides, C. Part III: Dynamic evolution of the particle size distribution in batch and continuous particulate processes: A Galerkin on finite elements approach. Chem. Eng. Sci. 2005, 60, 6998–7010. (10) McGraw, R. Description of aerosol dynamic by the quadrature method of moments. Aerosol Sci. Technol. 1997, 27, 255–265. 4001



ARTICLE

(11) Bekker, A. V.; Livk, I. An Implicit FEM Solution of a PBE Model of Gibbsite Crystallization with Constant and Nonlinear Kinetics. Ind. Eng. Chem. Res. 2011, DOI: 10.1021/ie1019326. (12) Li, T. S.; Livk, I.; Ilievski, D. Influence of the estimation procedure on the accuracy and precision of aluminium trihydroxide crystallisation kinetics from dynamic data. Ind. Eng. Chem. Res. 2001, 40, 5005–5013. (13) Safronov, V. S. Evolution of the Pre-Planetary Cloud and the Formation of the Earth and Planets [in Russian]; Nauka: Moscow, 1969. (Engl. transl. 1971 Israel program for scientific translations). (14) Bekker, A. V.; Livk, I. A PDE agglomeration model for Bayer precipitation. In Proceedings of the Chemeca 2009 Conference, Perth, Australia, Sept 27-30, 2009, [CD-ROM]; Engineers Australia: Australia, 2009; 10 pp. (15) Ilievski, D.; Hounslow, M. J. Agglomeration during precipitation: II. Mechanism deduction from tracer data. AIChE J. 1995, 41, 525– 535. (16) Misra, C.; White, E. T. Kinetics of precipitation of aluminium trihydroxide from seeded caustic aluminate solutions. Chem. Eng. Prog. Symp. Ser. 1971, 67 (110), 53–65. (17) Litster, J. D.; Smit, D. J.; Hounslow, M. J. Adjustable discretized population balance for growth and agglomeration. AIChE J. 1995, 41, 591–603. (18) Wynn, E. J. W. Improved accuracy and convergence of discretized population balance of Litster et al. AIChE J. 1996, 42, 2084– 2086. (19) Li, T. S.; Rohl, A. L.; Ilievski, D. Modeling non-stationary precipitation systems: sources of error and their propagation. Chem. Eng. Sci. 2000, 55, 6037–6047. (20) Li, T. S.; Livk, I.; Ilievski, D. Dynamic compartmental models of uniformly-mixed and inhomogeneously-mixed gibbsite crystallisers. Chem. Eng. Technol. 2003, 26 (3), 369–376. (21) Bekker, A. V.; Livk, I. 2-D Population balance modelling of gibbsite crystallization. In Proceedings of The 2007 AIChE Annual International Meeting, Salt Lake City, USA, Nov 4-9, 2007; AIChE: New York, 2007; 10 pp. (22) Alexopoulos, A. H.; Kiparissides, C. Part II: Dynamic evolution of the particle size distribution in particulate processes undergoing simultaneous particle nucleation, growth and aggregation. Chem. Eng. Sci. 2005, 60, 4157–4169. (23) Rosenberg, S. P.; Healy, S. J. A thermodynamic model for gibbsite solubility in Bayer liquor. 4th International Alumina Quality Workshop. Darwin, Australia; AQW Inc.: Australia, 1996; 10 pp. (24) Schibeci, M.; McShane, J.; Jones, W.; Vassallo, D.; Li, T.; Livk, I. Development of a laboratory semi-continuous gibbsite precipitator. 8th International Alumina Quality Workshop, Darwin, Australia, September 2008; AQW Inc.: Australia, 2008; 4 pp. (25) Bramley, A. S.; Hounslow, M. J.; Ryall, R. L. Aggregation during precipitation from solution: A method for extracting rates from experimental data. J. Colloid Interface Sci. 1996, 183, 155–165. (26) Li, T. S.; Livk, I.; Ilievski, D. Supersaturation and temperature dependency of gibbsite growth in laminar and turbulent flows. J. Cryst. Growth 2003, 258, 409–419.

4002


Comparison of FEM and DPB Numerical ... - ACS Publications

Recommend Documents