Mathematical analysis to optimize crystal growth

For the production of sodium sulfate, a brine is crystallized and crystals of glauber salt are generated by this process. The phase data related to the most common sodium sulfate minerals are as follows: mirabilite (Na2SO4 ⋅ 10H2O), tenardite (Na2SO4), glauberite (Na2SO4 ⋅ CaSO4), astrakanite (Na2SO4 ⋅ MgSO4 ⋅ 4H2O). The units commonly used to express the phases are moles of salt per 1000 moles of water. These latter units simplify the construction of the commonly employed four-sided Janecke phase diagrams. The cooling temperature or the speed with which the solution is cooled has an effect on the size and purity, as well as the amount of crystals produced. We seek to establish, through the population balance equations (PBE), which process variables can be modified to obtain a specific crystal size, as well as to validate the mathematical model that best predicts the amount of crystals precipitated as a function of temperature. The adjustment by least squares, cubic splines, pitzer equations and Lagrange interpolation is tested. The experimental results agree with the characteristics of the proposed models.


Introduction
Population balance equations (PBE) are widely used as modeling tools for particulate systems that estimate the dynamic evolution of particle size distribution (PSD) as a function of process operating conditions [1].The population balance is a continuity statement written in terms of the PSD and has consistently received attention since Smoluchowski (1917) [2] introduced the mathematical formalism over a century ago. A comprehensive overview of the mathematical issues involved, the numerical methods available, and possible developments for the future have been given by Ramkrishna (1985Ramkrishna ( , 2000 and Ramkrishna and Mahoney (2002) [3]. The population balance equations according to [4, Mullin 2001] express the relationship between crystal size L and the population density n (where G is the growth rate): We use this equation and we solve it by Gauss quadrature, with the objective of determining the moments μ 2 and μ 3 in order to know the growth rate of the crystals, which in turn are affected by τ (residence time), which is a process parameter that we can manipulate. On the other hand, we need to determine the optimal mathematical model to predict the precipitation of glauber salt as a function of the equilibrium temperature of the system (adiabatic crystallizers) since with the value (in g/L) of glauber crystals and the information obtained in the PSD we can estimate the weight of crystals in each fraction. With the mass of crystals in each fraction and the equation M i = f v ρL 3 , we can determine the population density, where f v is the shape factor of the particles [5].
With the aforementioned, it is very important to have a mathematical model whose predictive value is significant, since the crystallization potential (PCR g/L) is one of the variables that we use to obtain the growth rate in each of the crystallizers of the process.The most popular approaches for solving the PBE is the moment method, where internal coordinates are integrated out and governing equations are derived in terms of the desired moments. If the moment transform is applied to the PBE, it is possible to derive a continuity equation (i.e., the zero-order moment) and transport equations for mean properties [6]. We select the Gaussian quadrature method since it selects the evaluation points in an optimal way and not in an equally spaced way.
This equation contains 2n parameters to choose from. If the coefficients of a polynomial are considered parameters, the polynomial class of maximum degree 2n − 1 also contains 2n parameters [7]. For n = 2 in the interval [−1, 1] where the roots of the Legendre polynomials exist, the polynomial of maximum degree would be 2(2) − 1 = 3 such that: We use the roots for n = 2 of the Legendre polynomial that are x 1 = − √ 3/3 and x 2 = √ 3/3. The weights of the square c i are estimated by solving the linear system obtained by imposing the n conditions of accuracy, that is, the linear system: With this mathematical development, we seek to establish which process parameter can be used to manipulate the growth rate, in order to obtain crystals of desired size.
We can define based on this introduction the central aspects in which we work on the project. The optimal mathematical model is determined to predict the precipitation of the glauber salt as a function of the equilibrium temperature of the system (adiabatic crystallizers), which in turn determines which process parameter directly affects the growth rate through the population balance equations.
One aspect in the process control part that this work addresses is to integrate the mathematical model that predicts the amount of glauber precipitated in the control system of the production process, as well as to establish the mathematical model to determine the growth rate of the crystals.
Shows a schematic diagram on the most important aspects that are addressed in this work Fig. 1.

Prediction model for the precipitation of glauber salt
It is important to determine which of the mathematical models available and proposed in this study predict more accurately the precipitation of glauber salt as a function of temperature. To do this, glacier salt crystallization tests are performed at different temperatures; the temperatures selected to perform the test were 18, 15, 8, and 5 • C (international critical tables) [5]. I proposed four mathematical models of prediction to contrast them with the laboratory results; the selected models are least square method, cubic splines, change of coordinates (polar coordinates), and Pitzer equations.  In Fig. 2, we can see the tabulated data corresponding to the information contained in Table 1, on the values obtained in the crystallization tests and the values that the mathematical prediction models show.
The model chosen for the development of the tests is the one of cubic splines, the main reason for the choice of this model is that it provides an excellent fit to the tabulated points and its calculation is not excessively complex: for it, the nodes selected are 1 st node S 0 (t) (25,17.9 y 15.3); for the second node S 1 (t)(15. 3,8,5,0) is selected, where S is defined by a different cubic polynomial. Let S i be the cubic polynomial that represents S in the interval [t i , t i + 1], therefore: The polynomials S i−1 and S i interpolate the same value at point t i ; that is, it is fulfilled: Table 1 Crystallization potential values obtained from the mathematical models vs PCR values (g/L) obtained in the tests Mathematical prediction model   The mathematical configuration for the cubic splines in the nodes is: ⎛ The selection of the model was chosen based on the sum of the differences between the real value and that predicted by the mathematical model, according to the Table 2.
The crystallization tests are carried out at the indicated temperatures; we obtain the crystals at the laboratory level and we determine the mass of the glauber crystals at these reference temperatures (17.9,8,...,0) and compare them with the values obtained in the diagram of phases of Janecke. It is with respect to these data that the mathematical models mentioned above are proved.
The values obtained in the diagram (PCR g/L) are adjusted by means of a cubic polynomial, using the method of least squares, cubic splines, change of coordinates, and Pitzer equations, and are compared with laboratory values to determine which model reproduces more precisely what was obtained in the crystallization tests.

Mass balance
The change in dissolved solute between inlet and outlet of the vessel is matched by the gain in solid crystal mass. Thus:

Population balance
For a complete description of the crystal size distribution (CSD) in a continuously operated crystallizer, it is necessary to quantify the nucleation and growth processes and to apply the three conservation laws of mass, energy, and crystal population.
The crystal population density, n (number of crystals per unit size per unit volume of system), is defined by [4]: Where N is the number of crystals in the size range L per unit volume. The value of n depends of the value of L at which the interval dL is taken. The number of crystals in the size range L 1 to L 2 is thus given by: A population balance (input = output) in a system of volume V for a time interval t and size range L=L 2 − L 1 is: Where Q is the volumetric feed and discharge rate, G the crystal growth rate ( dL dt ), andn the average population density. As L → 0: Defining the liquor and crystal mean residence time τ = V Q and assuming that growth is independent of size ( L Law [4]), i.e., dG dL = 0, then: Which on integration gives: This equation is the fundamental relationship between crystal size L and the population density n characterizing the CSD.

Moments of distribution
Equation (4) is a number-size distribution relationship. Other distributions can be obtained in the following manner. From equations (4) and (5), the number of crystals, N, up to size L is given by: Average of the mass fractions μ 2 Average of the squared distance of possible values from the mean μ 3 Symmetry of the distribution (CSD)

Moment calculation by quadrature
In Gaussian quadrature, the points for evaluation are chosen in an optimal rather than an equally spaced way. The nodes x 1 , x 2 , ..., x n in the interval [a, b] and coefficients c 1 , c 2 , ..., c n , are chosen to minimize the expected error obtained in the approximation [7].
The coefficients c 1 , c 2 , ..., c n in the approximation formula are arbitrary, and the nodes x 1 , x 2 , ..., x n are restricted only by the fact that they must lie in [a,b], the interval of integration. This gives us 2n parameters to choose. If the coefficients of a polynomial are considered parameters, the class of polynomials of degree at most 2n − 1 also contains 2n parameters. This, then, is the largest class of polynomials for which it is reasonable to expect a formula to be exact. With the proper choice of the values and constants, exactness on this set can be obtained (Fig. 3). The quadrature approximation for the moment m of order k is [8]: Using the Gaussian quadrature for the interval of integration [−1, 1], thus: Where the coefficients A i and the points t i for i = 0, 1, 2, ..., n must be determined in such a way to get the best possible accuracy. For n = 2 in the interval [−1, 1] where the roots of the Legendre polynomials exist, the polynomial of maximum degree would be 2(2) − 1 = 3 such that: We use the roots for n = 2 of the Legendre polynomial that are x 1 = − √ 3/3 and x 2 = √ 3/3. The weights of the square A i are estimated by solving the linear system obtained by imposing the n conditions of accuracy, that is, the linear system: Solving the system we have to A 1 = A 2 = 1, then: Next, what we have to do is change the variable, to be able to work in the interval [−1, 1], i.e.:

Proposed solution and development
With the change of variables, we have: I k = Solution of the integral for the moments We use the solution proposed for moments and equation (3) to determine the growth rate with the moments μ 2 and μ 3 through the following: Where τ represents residence time.

Experimental control model for particle size in the crystallization process of glauber salt
Once the process parameter that allows me to modify the growth rate in the crystallization process is determined, a control method is proposed to give a residence time, to be able to obtain a desired size distribution (PSD), according to the following procedure: For a desired particle size distribution (PSD), the moments (μ k , k = 0, 1, ..., 3) obtained by solving the population balance equation and with the following relationship are used G = μ 3 μ 2 τ . We can obtain by the results of the crystallization tests at different residence times for each of the designated temperatures (18, 15, and 8 • C), the growth rate, and the residence times selected based on the feeding flexibilities of brine to the process. In addition, a very large residence time is selected in order to determine the effect on the particle size distribution, the residence time and the growth rate are plotted, and 6 points are obtained. Then, we use the Lagrange polynomial every 3 points to determine the growth rate at different residence times as shown in Fig. 4. P n = y 0 L n,0 (X) + y 1 L n, (X) + ...y n L n,n (X) Burden and Faires [7] Wehre L n,k is: And so, we can determine the growth rate, given a residence time contained in the interval the growth rate, as shown in the following table.
W e use the same criteria (interpolation polynomial) to nucleation rate (Fig. 5). By relating the D50 (particle size) to the growth rate of the crystals and adjusting the values with cubic splines, we can infer with the value of G interpolated by Larange in the previous graph (4) to obtain the desired crystal size as a function of the time of selected residence (Table 3).
With the six growth rate data obtained in the tests and with the Larange interpolation polynomial, we map the value of G for a considered residence time. This allows us to select the polynomial 1 or 2 of the cubic interpolation (6), and so determine the particle size for that considered residence time value. With the growth rate obtained given any residence time between the time interval in the crystallization tests, we can create logical conditions to select which cubic polynomial to use to determine the particle size, given a residence time, according to the following: Two interpolation polynomials are generated P 1,2 on the data of growth rate and residence time. Now, we select a residence time (τ ) value that is in the time interval considered in the tests and thus we obtain a growth rate G. With the growth rate obtained through interpolation, we evaluate in what part of the interval this value of G is understood, for this the interval is divided into nodes, every second node is evaluated if the value of G belongs to that interval. If this is the way, place the label of "True"; at the moment of appearing this label, we place the number 1 in the adjoining cell, and if this value appears, we map the coefficients of the cubic splines that relate the D 50 to the growth rate, as shown in the following table.
Below is the mathematical model of cubic splines for the speed of growth and the D 50 . This mathematical model is the one that is mapped in Table 4. We use cubic splines as an interpolation tool due to the distribution of points between the growth rate and the D 50 shown in Fig. 6. ⎛ We obtained similar values for the temperatures of 18 and 5 • C which indicate that the process parameter that we must manipulate is the residence time. This parameter is easy to manipulate since the crystallizers are of a constant volume, so by moving the flow of feeding we can manipulate the residence time at our convenience. Before being able to make these changes, it is necessary to determine with the current conditions what is the G value (growth rate) after determining at least another 2 points of G; these points should be at ends of the current condition (up and down), to be · · · · · · · · · · · · · · · · · · able to apply the interpolation polynomials and the coefficient selection criteria for a given τ to obtain a D 50 . This can be entered into the distributed process control and in this way determine the optimal D 50 values, as well as the mathematical model that determines the amount of glauber salt precipitated as a function of temperature (cubic splines).

Results
The tests carried out are the following: First, virgin brine is cooled to 18, 15, and 8 • C (values reported in international critical tables [5]) and at six different residence times: •  Fig. 4 is that what literation mentions is fulfilled [4], that is, at high residence times, the growth rate decreases. This can be explained according to the following: The pure condition of supersaturation or supercooling is not enough to cause a system to begin to crystallize. Before the crystals can develop, a number of small solid bodies, embryos, nuclei, or seeds, which act as crystallization centers, must exist in the solution. Nucleation can occur spontaneously or can be induced artificially. However, it is not always possible to decide if a system has nucleated on its own, or if it has done so through the influence of some external stimulus.As soon as a stable core has formed in an oversaturated system, they begin to grow to a visible size.
Diffusion theories assume that matter is deposited continuously on one side of the crystal at a rate proportional to the difference in concentration between the deposition point and the solution mass of solution, because a thermodynamic equilibrium is achieved since when the concentration between the solution and the crystals formed decreases, the rate of growth also decreases.
Below are the results for the test of 15 • C.

Conclusions
Through the development of this work, it is possible to establish a methodology through mathematical tools that allow the control of particle size distribution (PSD) as a function of a process variable. In this case, it is the residence time: this variable is controlled easily since the crystallizers are of constant volume, so the way to affect the residence time is by manipulating the volumetric flow of brine. It is also established which is the mathematical prediction model that best predicts crystallization as a function of temperature. It is important to mention that this tool was already programmed in the crystallization process and serves as a parameter to determine the growth rate.
The control logic to determine the PSD can be programmed in the distributed control system for the process at the industrial level and that this control is done online. In order to carry out the application, particle size measurement is needed in the process.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.