Dynamic simulation of industrial-scale gibbsite crystallization circuit

Population balance model is crucial for improving aluminum hydroxide massive crystallization and enhancing the control of precipitation trains. This paper presents the updated population balance model, which can be used for simulation of industrial-scale precipitation. Processes of birth-and-spread and particle breakage are considered integral parts of the precipitation process along with secondary nucleation, growth and agglomeration of particles. The conceptual difference of the proposed system of equations is its ability to reproduce the oscillatory process that occurs in precipitation circuits as a result of cyclic changes in the quality of the seed surface. It is demonstrated that self-oscillations can occur in the system without any external influence. The updated model is adjusted and verified using historical industrial data. The simulation of seed-recycle precipitation circuit showed an exact correspondence between the calculated dynamic pattern of changes in particle size distribution of aluminum hydroxide and the actual data.

and at certain points creates the risk of producing alumina, the granularity of which is not compliant with specifications.
The process data ( Fig.1, a) characterize the operation mode of precipitation areas at two alumina refineries, located in the Russian Federation: Urals (UAZ) and Bogoslovsk (BAZ) alumina refineries. The diagrams demonstrate the dynamics of fluctuations in the content of -45 µm fraction; regular changes of PSD can clearly be seen -with a period of around 2-3 months at BAZ and 1-2 months at UAZ. The amplitude of observed oscillations changes, their period is quite constant, and occurrence of the next extreme point seems inevitable.
Both refineries process boehmite-diasporic bauxites and use medium temperature digestion without sweetening. The caustic concentration in the pregnant liquor to precipitation amounts to 240-260 g/l (as Na2CO3); so, during precipitation A/С ratio in the liquor reduces from 0.55-0.57 to 0.28-0.31. Fine-seeded agglomeration is not very effective, therefore control of the precipitation process amounts to changing the temperature in head and tail precipitators, as well as altering the cones in product hydrocyclones (at UAZ).
One of the main objectives of such control method is to determine the best moment and level of involvement, so that control takes place in reversed phase to the natural flow of the process and its extent is adequate.
By using standard methods (microscope observations and laser analysis of particle sizes), it is possible to detect the coarsening or degradation after they have already started. In practice, attempts to take some preventive action prior to the appearance of instrumental signs of degradation or coarsening often resulted in a loss of control over the process. In this case, the peaks of -45 µm fraction content reached 1-3 % by the end of the coarsening phase, whereas in the degradation phase they run up to 60-70 %. In case of extensive coarsening, the yield of the area could decrease by 1.5 times. The changes in the content of certain particle grades, productivity, composition of the liquor, temperature and concentration profiles of precipitation at UAZ and BAZ refineries have been monitored for 15 years; process parameters are registered on a daily basis. It provided for the sufficient amount of data to study gibbsite crystallization under process conditions and to verify the simulation methods. Significant fluctuations, registered during observations, proved to be very useful for the present study.
A precipitation area at UAZ refinery was selected as the main object of research. At the refinery there are several precipitation trains of various capacity, but the general layout of their arrangement is identical (Fig.1, b). Precipitation of the liquor is conducted in a single stage. Batteries of hydrocyclones are used as classifiers. Only a part of the slurry flow is fed to classification in order to separate coarse fractions. The product hydrate, represented by hydrocyclone underflow, is filtered using vacuum filters and supplied to washing. The major share of the slurry flow from the last precipitator bypasses hydrocyclones and together with hydrocyclone overflow is fed to the station of seed filtration.
Population balance model. The fundamental equations used in the study are the equations of population and mass balance, which associate the precipitation level of aluminate liquor with the changes in the properties of aluminum hydroxide crystals [21]: where ni is the number of the ith grade particles; Li is the size of the ith grade particles, m; ∂n i /∂t the rate of change of particle number of certain size u over time in the slurry volume, 1/(m 4 s); G is the linear rate of crystal growth, m/s; Bi, Di are the rates of formation and disappearance of the ith grade particles as a result of their occurrence and transfer between grades, respectively, 1/(m 4 s); A is Al2O3 concentration in the liquor, g/l; M[Al2O3], M[Al(OH)3] are the molar masses of Al2O3 and Al(OH)3, respectively, kg/mol; А is the density of aluminum hydroxide, kg/m 3 . Models of growth and agglomeration. Equations (1)-(2) are solved simultaneously with expressions for crystal rates of linear growth and agglomeration, used in similar studies: • model of linear growth rate of the particles, independent of particle size [27]: where k0 is the pre-exponential factor of the linear growth rate constant; Eact is the activation energy of the reaction, J/mol; T is the temperature, K; Aeq is the equilibrium concentration of alumina in the liquor, g/l; z1 is the index of power at supersaturation; С is the caustic concentration of Na2CO3, g/l; • agglomeration model with a probability function, dependent on particle size, can be formulated for each particle grade as follows [11,15]: where v, u are the coordinates of particle in PSD; n is the number of particles in a volume; β(v,u) is the function, expressing the agglomeration probability of v-and u-sized particles (agglomeration kernel); k c is the proportionality factor in the agglomeration model; 4 is the correction to accommodate agitation conditions in the precipitator; G -particle linear growth rate. The first part in formula (4) describes the increment of particles in the analyzed v grade due to agglomeration of finer particles (from 0 to v), and the second part means the decrease of particles in the analyzed grade due to binding to agglomerates with other grades of particles.
Model of secondary nucleation. The rate of secondary nucleation or nucleation rate at the left boundary of the PSD (1/m 3 •s) is usually calculated at a constant rate of nuclei formation per unit area of phase contact [19]: where k n is the adjustable constant of nucleation rate, 1/(m 3 s); S sl is the phase interface area in the slurry volume, 1/m; z2 is the supersaturation exponent, in paper [19] z2 = 2.
Breakage model. The breakage of particles is usually associated with their attrition and cracking. Attrition results in insignificant changes in the size of the initial particles and formation of new nucleussized ones. In this sense, attrition is similar to contact nucleation so we will not analyze it separately and the term breakage means only an action of breaking of a particle to form new two similar particles, i.e. binary breakage.
In a general form the breakage function can be expressed as an agglomeration model for each particle grade [3]: where ϑ(t,u) is a function of breakage probability of uth grade particle at tth moment; the first part on the right-hand side of equation (7) represents the occurrence of uth grade particles due to breakage of coarser particles, and the addend corresponds to their breakage; n -number of particles in a volume. Simplification of equation (7) for binary breakage provides the following: Probability of particle breakage can be expressed by a power-law dependence [23]: where k b and p are adjustable constants used to consider the ratio of shear forces in the fluid and internal cohesive forces in particles;  is the shear rate, s -1 ; m 1 is the normalized first moment of crystal population in a volume, 1/m 2 . Hydrocyclone model. To complete the mass balance of precipitation along the seed classification contour, a hydrocyclone model proposed by Nageswararao was applied. The model was used to estimate classification size of the hydrocyclone: where K D0 is the correction factor for the classification size; Dc, Di, Do, Du are diameters of the cyclone, inlet, vortex finder and spigot, respectively, m; Lc is the length of cylindrical section of the cyclone, m; H is the hydraulic head, m; θ is the center angle of the conical section of the cyclone, deg;  is the hindered settling ratio, calculated by the following formula: where С v is the volumetric fraction of feed solids. The part of the liquid, discharged through the cone of the hydrocyclone (recovery of liquid to underflow), is calculated by the following formula: where K w0 is the correction factor. Fractional efficiency of particle capture was calculated using Whiten equation [2]: where , , ω * are adjustable constants. Oscillation character. There are different explanations for regular PSD fluctuations in seedrecycle precipitation circuits. Some authors point to the randomness of the process and fractal structure or roughness of the particles [28]; others associate the fluctuations with hydrate classifiers, which can cause PSD oscillations due to accumulation of aluminum hydroxide in the amounts, comparable to the solids amount in the precipitators [6]. The latter case, known in the theory of automatic control as a transport lag, under certain conditions can cause oscillations in the system, including non-damping ones. However, it does not explain the occurrence of a strong oscillatory process in crystallization circuits at UAZ, where seed classifiers are not used, and at BAZ, where aluminum hydroxide is not classified at all.
Under real process conditions, PSD fluctuations are also caused by perturbations of the flow rate, concentration, ratio and impurity content of the pregnant liquor [12], changes of daily ambient temperature, routine maintenance and shut-downs, and many other random factors. Nevertheless, no factor can explain the frequency, regularity and amplitude of the basic mode of oscillations shown in Figure 1, a.
A standard recommendation on nucleation control is monitoring and maintenance of a stable ratio of absolute supersaturation to full seed surface С/Ssl [26]. Performed statistical analysis of industrial data from UAZ and BAZ refineries showed that the correlation between the rate of PSD changes and the ratio of supersaturation to full seed surface С/Ssl is very weak (determination factor R 2 is less than 0.1). Therefore, the recommendation is not applicable for these refineries.
As a working hypothesis on the character of oscillations, an assumption was adopted in the present study that the driving source of the longest wave of PSD fluctuation has natural causes and is associated with seed surface properties. This assumption complies with the generally accepted opinion that the moment, when the nuclei start actively passing into the liquor, is associated with a decrease in number density of linear growth centers on the surface of the crystal [7].
Birth-and-spread model. According to crystallization theory, the processes of growth and nucleation have a common beginning, which involves molecule adsorption on the surface of a particle with partial loss of the solvation sphere (Fig.2, a). Initial attachment of the molecule occurs on a flat, energetically least favorable site of the particle (Fig.2, b). Then the absorbed molecule migrates along the surface and overcomes energy barriers 3 and 5 to reach more favorable sites for fixation 4 and 6 [7].
In case route 3 is unavailable, ripening of the nucleus and addition of new molecules will occur on the flat site; in this case the bond between this group of molecules and the base one will be the weakest, and subsequent separation of the nucleus by external mechanical force is most likely to happen. Step site Flat site Taking into consideration the indicated mechanism, the adsorption process of embryos on the seed surface does not depend on surface quality and is determined only by its area, diffusion factor, diffusion layer thickness and liquor supersaturation. Thickness of the diffusion layer can be considered conditional-constant at a fixed agitation rate. According to Stokes-Einstein formula, the diffusion factor depends on temperature, but there is no need to consider its impact separately, as in the first approximation this impact can be taken into account when tuning the index of power during supersaturation. Hence, to calculate embryo adsorption rate on the seed surface, equation (6), proposed for secondary nucleation, can be used without alterations.
Energy consumption, associated with the embryo navigating routes 3 and 5 ( Fig.2, a, b), characterizes the probability of its transfer from seed surface to the liquor. If these routes are shortened, mainly due to agglomerate formation, increasing surface fractality and leading to generation of numerous growing sites on the crystal, the probability of nuclei fixation on the seed surface increases, while the inverse probability of nuclei passing into the liquor decreases.
As the particle grows and its shape improves, these routes get longer, embryo fixation in the growing points of the parent particle deteriorates, and the embryo has more chances to form its own growing point and to pass into the liquor as a nucleus.
It is proposed to describe birth-and-spread as a two-stage process, during which an embryo is formed on the surface and then grows together with the particle, if suitable vacant sites are available near the fixation site (Fig.2, c). The process rate is proportional to the product of two probabilitiesthe probability of embryo absorption on the seed surface and the probability of its firm bonding with the surface. Therefore, the expression for birth-and-spread rate will be written as follows: (14) or where k 1 is a dimensionless coefficient, considering other factors; Nv/Ne = nv is the number of available vacant sites for one embryo; μ i = nv/S f is the number density of available vacant sites per unit surface of the particle, 1/m 2 ; ∑ S f.i μ i i is the average surface vacancy over total particle population. The number density of available vacant sites μ i is an expression of seed surface quality. It decreases as a result of prolonged growth of the particle and increases during agglomeration. Consequently, for the population of particles with prevailing agglomerates it will be higher than for the population with a small fraction of agglomerates. In this way, the fraction of agglomerates can be related to the number density of vacant sites. In order to avoid using absolute values of the number of vacant sites, a transition to relative values is required: where μ а and μ g are number densities of available vacant sites for agglomerates and growing particles, respectively, 1/m 2 ; f a and f g are the fractions of agglomerates and growing particles, respectively; μ а /μ g =  is the value expressing how much the probability of an embryo binding with the agglomerate surface exceeds the same probability for the growing particle.
Assuming μ g to be constant, expression (15) for birth-and-spread rate with due account of formula (16) can be rewritten as follows: where ks = k1g is an adjustable constant. The fraction of agglomerates fa is a relative notion, as the same particle can be a newly formed agglomerate of two particles or an agglomerate, which for a long time has been growing as a single particle. For the sake of consistency, the authors suggest that particles are divided into agglomerates and non-agglomerates in the moment of their transfer between size grades. At each iteration, when new particles emerge in the grade, it will enable to register their type and calculate a new fraction of agglomerates in the grade. If the particle in the analyzed grade is formed as a result of nucleation, linear growth or breakage, it is classified as a common particle; if its formation is associated with agglomeration process, it is treated as an agglomerate.
Discussion. Updated model of secondary nucleation. Taking into consideration the competition between secondary nucleation and birth-and-spread processes, with due account of formulas (6) and (17), the expression for the rate of nuclei passing into the liquor can be formulated as follows: or where B 0 is the rate of nuclei transition to the liquor, 1/(m 3 s); B M is the overall rate of nucleation on the particle surface, calculated using formula (6), 1/(m 3 s); B S is the rate of birth-and-spread process, 1/(m 3 s); B M /B S is the ratio between the rates of secondary nucleation on the particle surface and birthand-spread process. k n , k s and  parameters in formula (19) are adjusted in a dynamic mode in the process of model tuning by simulating actual PSD oscillations.
Updated breakage model. Equation (8) adopted the binary pattern of particle breakage. In order to consider the possibility of multiple sequential breakages under such conditions, the breakage probability function for particles of different fractions needs to be introduced. For this purpose, a Gaussian function was used in the study, the peak of which corresponds to the size of particles that break the most: where F(u) is the breakage intensity of the uth grade particle, 1/(m 3 s); L * is the size of the particles that break the most, m. The presence of multiple impurities hinders layer-by-layer linear growth of the particles, which causes formation of irregular particles and dendrites, poorly connected to the basis and prone to shearing. Such particles form during the period of active agglomeration, especially if it occurs at low temperatures. Similar to the birth-and-spread model, in order to consider the intensity of particle breakage in each grade, the fraction of agglomerates ( , ) can be used. After supplementing the equation of particle breakage probability (9) with new parameters, the following equation is obtained: Thus, the proposed updated model of particle breakage consists of equations (8), (20) and (21). Oscillatory process features. The introduced equations, comprising agglomeration model (4)-(5), updated nucleation model (19) and updated breakage model (8), (20) and (21), can be simplified and combined in a separate subsystem of ordinary differential equations: where Ṅg .a = k ga N g (t)N a (t) + 2k gg N g 2 (t); (24) Ṅa .a = k gg N g 2 (t) -k aa N a 2 (t); where Ṅg .n , Ṅg .a , Ṅg .b are terms, indicating the change rate of growing particle number due to nucleation, agglomeration and breakage, respectively; Ṅg .b , Ṅa .a , Ṅa .b are terms, indicating the change rate of agglomerate number due to nucleation, agglomeration and breakage, respectively; N g (t), N a (t) are the numbers of growing particles and agglomerates, respectively; f a = N a (t)/ N g (t) + N a (t) is the fraction of agglomerates in the population; [N g (t) + N a (t)] is a simplified assumption for full surface of the particles Ssl in formula (19); k n1 , k n2 , , k ga , k gg , k aa , k b are constant factors.
The system of equations (22)-(27) is solved under initial conditions dN g (0) = n g and dN (0) = n a . Mass sources and sinks were excluded from the system, liquor properties were assumed constant, particle growth was not taken into account. Under such conditions, the system can be analyzed for self-oscillations only in respect to seed surface properties.
The structure of equation system (22)-(27) is similar to a well-established Lotka-Volterra model of predator-prey relationships [17], only in the course of precipitation the competition takes place between growing particles and agglomerates as a result of secondary nucleation, birth-and-spread, agglomeration and breakage processes.
Numerical simulation. Numerical simulation is performed using PrecipExpert software, integrated with the obtained model [1]. The software enables to create flow diagrams of industrial-scale precipitation from separate units, locate them on the worksheet and connect in the required order. The system of population balance equations is solved using discrete method [10,14].
Numerical simulation were performed for the precipitation flow diagram shown in Fig.4. It comprises 16 precipitators Рис.3. Self-oscillatory process in a closed loop system during precipitation P6L1_1-P6L1_16 with the volume of 3500 m 3 , a battery of hydrocyclones CN_6, a seed filter VFz_6, a product separation filter VF1p_6, a first stage washing filter VF2p_6 and agitators. Pregnant liquor and seed, inherent to the current precipitation train, are fed to the mixer TKh_6_1, the seed from another precipitation train is fed to this mixer. Additionally, the flow of liquor and finely-dispersed aluminum hydroxide CleanLiq_P6, generated in the process of chemical cleaning of precipitators, is considered.
Different periods of precipitation train operation were selected for training and prediction. Initial data were taken from the historical database of the refinery and additional reports, which were added into the model in a dynamic mode during the solution process. The initial data were updated once a day, daily average values were used. Preparation of the initial data involved only removing outliers; data smoothing was not applied. The following parameters were loaded to the model as initial data: flow rate and composition of the pregnant liquor, solid concentration in the first precipitator, temperature in the precipitators, number of hydrocyclones in operation, spigot diameters of hydrocyclones, flow of additional seed from another precipitation train, volumes and compositions of liquors, supplied after cleaning of precipitators.
Data on aluminum hydroxide PSD were introduced to dynamic calculations only once, at the beginning. Then followed a short period of model self-tuning. The rest of the time, the PSD of aluminum hydroxide in each precipitator and the station was calculated by the model along with liquor compositions in all intermediate and outlet flows.
The system of population balance equations was solved using a mesh, comprising 30 grades: Flow exchange period between the units (external transportation problem) amounted to 10 minutes; the population balance model was solved at a time step of 150 s. Execution time of one dynamic calculation over the period of 1.5 years was approximately 8 hours using a CPU with 133 GFlops.
The tuning of the model was performed in two stages. At the first stage, factors of growth, agglomeration and breakage models were selected for average values of all input parameters under steady-state conditions. In the growth model (3), parameter was adjusted in such a way that the estimated ratio of the pregnant liquor corresponded to the average value for the monitored period, when the process reached a stable state. In a similar manner, the average content of -1 µm, -5 µm, -20 µm, -45 µm particles were adjusted by selecting factors in formula (5), in formula (20), and p in formula (21). The contents of the coarsest grades +100 µm and +125 µm were adjusted by changing hydrocyclone parameters and the curve of fractional particle capture efficiency (10)-(13).   The second stage of model tuning comprised the selection of k n , k s and μ factors in formula (19) in a dynamic mode. Several calculations were performed in parallel for different factor values. Optimum values were determined by means of comparing calculated curves of particle size changes with actual process data. Least square method was used for quantitative estimation of the simulation quality. As a result of model tuning, the following set of factors for the model was obtained: kn = 330; Eact = 6.3210 5 ; z1 = 4; ks = 0.61; kb = 0.95; p = 1;  = 3.33;  = 2.5; kn = 110 14 ; z2 = 2; ks = 12.
Quality tuning was also performed for the initial model, in which nucleation was described using formula (6), and breakage was described using formulas (8) and (9).
Results. The initial population balance model for precipitation and the updated model with suggested modifications -the equation (6) in nucleation model was replaced with (19), the equation (9) in the breakage model was replaced with (21) and supplemented with equation (20) -were solved under the same initial conditions and environmental variables. Duration of the simulation period amounted to 17 months. Over this period, three complete cycles of PSD fluctuations were registered.
The dynamics of changes in two most important particle grades -20 and -45 μm, calculated using the models and compared to the process data, is presented in Fig.5, a, b. Over the analyzed period, along with three maximum and four minimum peaks, the content of the main fractions demonstrated a large amount of minor fluctuations, which can be attributed to changes in the external factors. These factors produced a significant impact on both calculations. The prediction, obtained by the initial model, was not satisfactory. Calculated and actual waves were out of phase, and the amplitude of predicted fluctuations was approximately two times smaller than the actual one. The updated model reproduced the dynamics of PSD fluctuations much more accurately. For practical purposes it is important that the updated model can accurately predict the occurrence of PSD peaks and their amplitude. The Table specifies  Deviation in the peak date of the coarsening phase is three days, of the degradation phase -7 days.
Significant discrepancies in the content of -45 µm fraction between new model calculations and actual data in January and September 2015 are explained by the fact that prior to these two periods the system was subject to a strong external impact -within 1-2 days seed concentration in the slurry decreased sharply from 820-850 to 450-500 g/l. After that, fine seed from another precipitation line was added to the process; under the conditions of the problem, it did not contain any agglomerates. The model reacted to these actions by a consistent decrease in particle breakage rate, increase in agglomerate content, growth of secondary nucleation and, consequently, a reduction in the number of new nuclei passing into the liquor (Fig.5, d).
After failed results of simulation in January 2015, the new model corrected itself by demonstrating the next maximum peak of -45 µm fraction content in February and the minimum peak in April 2015.
The best agreement between the calculation results and the observations was achieved for A/C ratio in the spent liquor (Fig.5, c). Strong noise in the diagram is explained by fluctuations of composition and flow rate of the aluminate liquor and addition of seed from other areas, whereas the apparent wave is associated with changes in specific surface area of the intrinsic seed.
According to simulation results, the majority of agglomerates occur within the size ranges of 0.5-1.0 and 4-16 µm (Fig.6, a). The presence of the first peak aligns with the notion that 0-2 µm particles demonstrate an intensified growth compared to the remaining population [7]. Constant presence of the second peak is registered by means of microscopic investigation, which proved that particles of approximately 10 µm size are characterized by strong fractality (Fig.6, b).
According to simulation results, the agglomerate fraction in particle grades is not constant and subject to fluctuations just like other properties of particle population. In the periods of seed coarsening the agglomerates get coarser too, which causes the first peak in their content to disappear. In the periods of seed degradation, the agglomerate fraction in finer particle grades increases (Fig.6, a).
As a result, the updated model after performed tuning demonstrates a good predictive power for the content of coarse particles -the correlation ratio is within the range of 0.78-0.85, and mean square deviation, e.g. for -45 µm fraction, amounts to 6 %. As the particles get finer, predictive accuracy decreases. To control the precipitation process, it is critically important that the use of the updated model enabled to achieve adequate predictive accuracy for the coordinates of the next inflection point in the fluctuation curve of the -45 µm fraction content (see Table). Conclusions.
1. Basing on the performed study, a conclusion can be drawn that the birth-and-spread process is one of the forms of particle population development during precipitation of pregnant liquor at low temperatures. It is governed by seed activity, accelerating simultaneously with agglomeration and slowing down with the reduction of agglomerate fraction in the population. Birth-and-spread serves as an inhibitory process, which hinders the ingress of secondary nuclei into the liquor.
2. The competition between the processes of secondary nucleation, agglomeration and breakage of particles gives them properties of an oscillatory system. These oscillations are stable and under certain conditions can become non-damping. To determine fixed points, oscillation phases and the rate of exciting processes, a new property of the population can be introduced -the fraction of particles represented by agglomerates in different size grades.
4. The absence of parameters characterizing seed surface quality in the equations (6) and (9), which express the intensity of nucleation and particle breakage, limits their application to crystallization circuits with fine seed agglomeration at high temperatures. Controlled agglomeration enables to supply active seed to the second stage of precipitation and to control the population of agglomerates. This significantly decreases the amplitude and oscillation period for PSD and liquor productivity. In single-stage precipitation circuits at low temperatures it is impossible to maintain constant activity of the seed surface, and that must be considered in mathematical modeling.
5. The proposed updated models of nucleation, birth-and-spread and breakage processes cannot fully explain all PSD changes, observed during precipitation, as well as available historical data on the process do not reveal all the variability of external conditions. Nevertheless, obtained calculation results demonstrate that the model quite accurately reproduces the phase, amplitude and period of the basic mode of oscillations.