CET Volume 86 DOI: 10.3303/CET2186146 Paper Received: 12 August 2020; Revised: 28 February 2021; Accepted: 9 May 2021 Please cite this article as: Elmisaoui S., Latifi A.M., Khamar L., Salouhi M., 2021, Shrinking Core Approach in the Modelling and Simulation of Phosphate Ore Acidulation, Chemical Engineering Transactions, 86, 871-876 DOI:10.3303/CET2186146 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 Shrinking Core Approach in the Modelling and Simulation of Phosphate Ore Acidulation Sanae Elmisaouia,b,c, Abderrazak M. Latifia,b,*, Lhachmi Khamara,d, Mohamed Salouhic a Mohammed VI Polytechnic University, Benguerir, Morocco. bLaboratoire Réactions et Génie des Procédés, CNRS-ENSIC, Université de Lorraine, Nancy, France. cL3GIE, EMI, Mohammed V University, Rabat, Morocco. d LIPIM, National School of Applied Sciences in Khouribga, Sultan Moulay Slimane University, Beni Mellal, Morocco. Abderrazak.latifi@univ-lorraine.fr This paper deals with the modeling and simulation of the acidulation (dissolution) of phosphate ore particles in a dilute phosphoric acid solution in a batch stirred tank reactor. A shrinking core model with elimination of the products is used to describe the reaction and mass transfer phenomena involved in the three phases considered, i.e., liquid bulk, liquid film surrounding the particles and solid phase. The model is based on mass balance equations in the three phases and the corresponding equations are ODEs in the liquid bulk and in the solid phase, and PDEs in the liquid film. An estimability analysis method based on global sensitivities is used to determine the most estimable unknown parameters involved in the model equations from the available experimental data. The experiments consist of measurements over time of the conversion rate of the ore particles obtained in operating conditions similar to those used in an industrial phosphoric acid plant. The values of the least estimable parameters are fixed from the literature, while gProms environment is used to implement and solve the equations and to identify the most estimable parameters. The results obtained show a good agreement between the model predictions and the measurements with coherent values of the identified parameters. 1. Introduction In the phosphate industry, phosphate ore is often mined from open-pit mines and then beneficiated in a series of flotation units. The ore thus upgraded is then sent to a (wet) phosphoric acid process to produce phosphoric acid which is the key component of the phosphate industry where it is used in the production of fertilizers, in the food industry and even in the pharmaceutical industry (Becker, 1983). The manufacturing process which is widely used in phosphate industry consists mainly of a digestion tank, a filtration unit, and a concentration unit. However, despite its wide use, its functioning still poses several problems due the lack of knowledge on many complex phenomena involved in the process. The understanding of these phenomena and their fine modelling along with accurate experimental measurements of the relevant variables are some of the key issues that should be addressed to optimally design and operate the phosphoric acid processes. Acidulation (or dissolution) takes place in the digestion tank and consists in reacting the phosphate ore with a mixture of concentrated sulfuric acid and recycled phosphoric acid to produce phosphoric acid and gypsum. The digestion tank which is the heart of the process is therefore critical for the performance of the downstream units, i.e., filtration and concentration units. The control of its operation is therefore of utmost importance for the entire phosphoric acid process. The objective of this paper is to analyze the acidulation mechanism of the phosphate ore particles in the digestion tank. The acidulation is considered as a set of non-catalytic fluid/solid chemical reactions with consumption of the solid. Its modeling is therefore based on the approach of shrinking core model with elimination of the products. Furthermore, the model involves several phenomena including mass transfer and transport of reagents (products) to (from) the particle surface, and reactions at the surface of the solid. The model equations describe the mass balances within the liquid bulk, the liquid film surrounding the particles, 871 and the particles, and involve several unknown parameters. A global estimability analysis is then carried out to determine the parameters which are estimable from the available measurements. The most estimable parameters are then identified from the measurements of the phosphate ore conversion rate, and the values of the less estimable parameters are fixed from the literature. gProms environment will be used to implement and solve the model equations and for parameter identification as well. 2. Model formulation Tri-calcium phosphate (TCP) is one of the most interesting elements to be recovered from the phosphate ore. The objective of the acidulation is to dissolve the TCP in a dilute phosphoric acid solution before its attack with a concentrated sulfuric acid solution to produce concentrated phosphoric acid and gypsum. Two main reactions take place in the digestion tank during the dissolution process (Elmore and Farr, 1940) as: Ca (PO ) + 4H PO → 3Ca(H PO ) (1) Ca (PO ) + H PO → 4CaHPO (2) The conversion rate of those reactions depends on different operating parameters, including temperature, phosphoric acid concentration, porosity and size distribution of phosphate ore particles, residence time, solid rate and hydrodynamics (Theodosiadis et al., 2007). In this work, we will focus on the most important one which transforms TCP into mono-calcium phosphate (MCP) (Eq.1). 2.1 Dissolution mechanism The attack of the phosphate ore particles by a dilute phosphoric acid is a non-catalytic liquid/solid reaction with a complex kinetic mechanism. We assume that three phases are involved in the mechanism, i.e., liquid bulk, liquid film surrounding the particles, and solid phase (Figure 1). Figure 1: Schematic illustration of the SCM in the dissolution mechanism The dissolution mechanism considered assumes that the phosphoric acid (H+ ions) transfers from the liquid bulk to the liquid film, then diffuses through the film towards the solid, adsorbs to the phosphate ore particles, and finally reacts with the TCP reagent at the solid surface. The MCP product takes the opposite path where it first desorbs from the solid surface, then diffuses through the liquid film towards the interface, and finally transfers to the liquid bulk. It is noteworthy that all the mass transport/transfer phenomena and chemical reaction take place simultaneously in the digestion tank. 2.2 Rate limiting step The digestion mechanism described above is usually controlled by a rate limiting step. The latter can be: • the transfer of H+ ions and their diffusion from the solution bulk to the liquid-solid interface (Ben Brahim et al., 1999; Jamialahmadi et al., 1998), • the chemical reaction of the phosphate ore particles with the H+ ions, • and the transfer of products from the liquid-solid interface to the bulk solution (Van der Sluis, 1987). However, the rate limiting step is still controversial and not yet known precisely despite the large number of research contributions devoted to that issue which deserves to be further investigated. 872 3. Phosphate ore dissolution modelling The dissolution model developed is based on the following assumptions: • Particles are well dispersed in the liquid phase, are spherical of the same diameter, and shrink uniformly during the dissolution process, • Digestion tank is perfectly mixed, • Reaction is irreversible and only takes place at the surface of the particles, • Adsorption and desorption steps are assumed to be very fast and are not considered. Moreover, the shrinking core model (SCM) with elimination of the products used is adapted from the model developed in (Salmi et al., 2017) with two main major differences. The first one assumes that the solid reagent is not soluble in the liquid, which is the case of the Moroccan phosphate ore, and the second one considers that the surface of the particles is not saturated with the solid reagent. The model equations are based on the mass balances in the liquid film, in the liquid bulk and in the solid phase. They are presented below. - Considering the realistic pseudo steady-state in the liquid film, the mass balance can be written as: 1r ∂∂r D r ∂C∂r = 0 ; i = H PO , MCP (3) with the following boundary conditions: •at the solid interface: − D = ℜ (4) where the kinetic rate equation is assumed to be of first order with respect to phosphoric acid at the particle surface and defined as: ℜ = ν k C , kr is the rate constant and is the stoichiometric coefficient of component i in the reaction (1). •at the liquid film / liquid bulk interface: − D = k C | − C′ , (5) The equations (3), (4) and (5) allowed us to express the components’ concentrations at the surface of the phosphate ore particles as follows: C = (6) C C + 3 C (7) - The mass balance in the liquid bulk is expressed as: = −a D (8) -and the mass balance in the solid phase is given by: ( ) = ( C′ − C | ) (9) The initial conditions are given in Table 1. Table 1: Initial conditions used in the simulation Phase H PO Ca (PO ) MCP Liquid film C (mol/m3) 0 — 0 Bulk solution C (mol/m3) 180 — 0 Solid n ( ) (mol) — 0.0322 — The liquid film thickness is calculated using a mass transfer correlation as (Salmi et al., 2017): 873 δ = R 1 + α RR D (10) and ′ are the concentrations of component i in the liquid film and in the liquid bulk, respectively, and are the diffusion coefficient and the transfer coefficient of component i, respectively, n ( ) is the mole number of TCP in the solid phase, is the initial radius of the particle, is a hydrodynamic parameter, and is the particle specific area, i.e., particle area to liquid volume ratio at the beginning of the reaction. The assumption of spherical particles leads to the following relation between their specific area and radius: = . The dynamic model is described by a set of partial differential equations (PDEs) in the liquid film and ordinary differential equations (ODEs) in the solid phase and in the liquid bulk. The model equations are implemented and solved within gProms software which provides an environment for modeling, simulation, and parameter identification of dynamic systems (Process Systems Enterprise, 1997-2020). Furthermore, the unknown parameters are identified from the data available in the literature and obtained in operating conditions similar to those used in an industrial phosphoric acid plant. 4. Results and discussions Figure 2: Example of concentration profiles of phosphoric acid (a) and MCP product (b) in the liquid bulk, solid conversion rate (c), and liquid film thickness (d), obtained at different hydrodynamic conditions. An example of the simulation results is presented in Figure 2. It shows the concentration profiles of acid and MCP in the liquid bulk (Figures 2a and 2b), the conversion rate of the solid phase (Figure 2c), and the liquid film thickness (Figure 2d), obtained at different values of the hydrodynamic parameter α which is a function of the stirring speed. For a given value of α, the concentration of acid decreases over time in the bulk solution (Figure 2a). At the same time, the MCP produced diffuses from the surface of the particles to the bulk of the liquid solution, and its concentration increases with time (Figures 2b). Finally, Figures 2c and 2d show that the TCP conversion rate increases, and the liquid film thins as the phosphate particles shrink, thus showing that the phenomenon of digestion takes place. Furthermore, as expected, when the stirring speed increases, and consequently α, the phosphoric acid concentration as well as the liquid film thickness decrease, while the MCP concentration and the solid conversion rate increase. (a) (b) (c) (d) 874 On the other hand, as can be noticed, the model involves several unknown parameters (kr, α , D ,D , k and k ) that should be identified from the available experimental measurements by means of a parameter optimization method. According to the film theory, the transfer coefficient can be expressed as a function of the diffusion coefficient as: k = D δ⁄ , where is the liquid film thickness. The number of unknown parameters is therefore reduced to four, i.e., kr, α , D and D . In this work, the available experiments consist of measurements of the TCP conversion rate over time in a stirred batch reactor where the digestion of the phosphate ore particles by a dilute phosphoric acid solution takes place (Ben Brahim et al., 1999). Thus, using our recently developed estimability analysis method based on global sensitivities (Bouchkira et al, 2021) to rank the unknown parameters from the most to the least estimable, we were able to assess whether the experimental data available contain the information necessary to identify all unknown parameters of the acidulation model. 5. Parameter identification The use of the aforementioned estimability analysis method leads to the following order of estimability of the unknown parameters: kr > > > . Then, a gradient-based optimization method was used to identify the most estimable parameters, i.e., kr and , while and are set to values taken from the literature (Ben Brahim et al., 1999). The values of the identified parameters and their corresponding 95% confidence intervals (CI) are summarized in Table 2. Table 2: Identified parameters values Parameter Value 95% CI Unit 8.68x10-5 ±1.05x10-5 m/s 4.98x10-2 ±1.39x10-2 ----- The tight confidence intervals show that the parameters are accurately determined. Moreover, Figure 3 exhibits a good agreement between the model predictions of the conversion rate of phosphate ore particles and the measurements carried out at a temperature set at 75°C. Figure 3: Comparison of simulated and measured values of solid conversion rate Figure 4: Evolution of particle radius, film thickness and transfer coefficients with time As the dissolution proceeds, the radius of the particles as well as the thickness of the film decrease, while the mass transfer coefficients of acid and MCP increase over time (Figure 4). They all reach a plateau after only 400 s of dissolution, as does the conversion rate (Figure 3). On the other hand, the optimal value of the hydrodynamic parameter α is exploited to determine the stirring speed used to carry out the experiments presented in Figure 3. Indeed, the mass transfer correlation used to express the film thickness Eq. (10) enables us to express the hydrodynamic parameter as: 875 α = (11) where is the energy dissipated in the reactor given as: = (12) and are the kinematic viscosity of phosphoric acid and the density of the reaction mixture respectively, is the stirring speed, is the power number, D is the agitator diameter, and is the initial mass of the particles. The parameters used to perform the simulations are shown in Table 3. Table 3: Model parameters Parameter Value Unit 9.6x10-7 m2/s 986 kg.m-3 3x10-1 --- D 7x10-2 m 1.5x10-4 m 10-2 kg The stirring speed deduced from Eq. (12) is about 796 rpm which is almost the same as the value of 810 rpm used in the experiments (Ben Brahim et al, 1999), thus showing that the model developed can accurately describe the dissolution mechanism of the phosphate ore in a dilute phosphoric acid solution. 6. Conclusions A first-principles model is developed and allowed us to understand some of the complex mechanisms of the dissolution of phosphate ore particles in a dilute phosphoric acid solution. However, many issues are still to be addressed and developed to improve the prcision of the model predictions. They include the phenomena involved at the surface of the particles, i.e. adsorption of the reagents and desorption of the products, and the kinetics of the reactions involved. Furthermore, the rate limiting step of the dissolution phenomenon as well as the chemical equilibrium deserve to be further investigated. Moreover, additional digestion experiments should be carried out to accurately calibrate and validate the model prior to its use in subsequent optimization of the digestion tank performances. References Becker P., 1983, Phosphates and phosphoric acid: Raw materials, technology and economics of the wet process, Fertilizer Science and Technology Series. New York. Ben Brahim F., A. Mgaidi, D. Oulahna and M. El Maaoui, 1999, Kinetics of leaching of tunisian phosphate ore particles in dilute phosphoric acid solutions, The Canadian Journal of Chemical Engineering, 77, 136-142. Bouchkira I., A. M. Latifi, L. Khamar and S. Benjelloun, 2021, Global sensitivity based estimability analysis for the parameter identification of Pitzer’s thermodynamic model, Reliability Engineering and System Safety,207, 107263. Elmore K. L. and T. D. Farr, 1940, Equilibrium in the system calcium oxide-phosphorus pentoxide-water, Industial Engineering Chemistry, 32, 580-586. Jamialahmadi M., S.H. Emam and H. Miiller-Steinhagen,1998, Dissolution of Phosphate Rock by Mixtures of Sulfuric and Phosphoric Acid, Developments in Chemical Engineering and Mineral Processing, 273-293. Pitzer K. S., 2017, Activity Coefficients in Electrolyte Solutions, 2nd edition, CRC Press. Salmi T., V. Russo, C. Carletti, T. Kilpiö, R. Tesser, D. Murzin, T. Westerlund and H. Grénman, 2017, Application of film theory on the reactions of solid particles with liquids: Shrinking particles with changing liquid films, Chemical Engineering Science, 160, 161-170. van der Sluis S., 1987, A Clean technology: Phosphoric acid process. Delft University Press. Theodosiadis K., A. I. Papadopoulos, P. Seferlis, 2017, Modeling of the reactor-crystallizer unit in a phosphoric acid production plant, Chemical Engineering Transactions, 12, 237-242. 876