domain: spi info.info-ia phys

text:
Safae El Misaoui

M Abderrahman Nounah

M My

Saddik Kadiri

M Khalid

Achoubir Pr

M M'barek Bakkas

Jean-Michel Ghidaglia

M Radouan

Boukharfane Pr At

M Saad Benjelloun

M Lhachmi Khamar

M Mohamed Khamar

CFD modeling and simulation of complex flows in phosphate fertilizer manufacturing process: Preneutralization and granulation units By

Keywords: Ammonium Phosphate Fertilizer, Preneutralization, Granulation, Multiphase Reacting Flows, Modeling, Simulation and Optimization viii

Mathematical modelling and numerical simulations are widely used in many chemical industries, such as petrochemical industry. Its used for operator training, design, and process optimization. However, there is a lack of rigorous numerical modeling and simulations in the phosphate fertilizer industry. There exist many challenges in the production systems of phosphate fertilizers including; multiple phases in the system such as liquids, solids and gas, particles with different sizes and distributions, and variation in the physical properties including rheology and thermodynamics properties. Despite the available resources, limited studies have been focused on the understanding of the hydrodynamics within the components of the manufacturing process, to enhance the production yield and to adapt the manufacturing process conditions to the quality of raw materials. Computational Fluid Dynamics (CFD) represents a valuable tool to validate phenomenological models and allows an exhaustive understanding of process sensitivity to operating conditions, and geometrical details effects on the flow hydrodynamics. Furthermore, CFD permits a better understanding of physical phenomena, as exchange mechanisms occurring in the studied processes, which is not possible on the industrial plant or by laboratory scale experiments that require many considerations for the scale up. The aim of this thesis is to i) investigate the fluid dynamics of the master pieces in Di Ammonium Phosphate fertilizer process, which are the preneutralization and the granulation unit operations, using commercial simulation software Ansys Fluent, ii) to identify the effects of the complex phenomena in the formation of heterogeneities of the multiphase flows, and iii) to evaluate the effect of the operating conditions on the flows hydrodynamics. On the light of the obtained results, we aim to optimize the design and control of the preneutralization and granulation units and propose new recommendations for a better process performance. This research work is composed of two different studies; on the one hand, the three-dimensional modeling and simulation of the preneutralizer reactor, based on the Euler-Euler multiphase approach, and the Multiple Reference Frame (MRF) is applied to consider the rotational movement of the agitator. The tailored gas-liquid CFD model accounts for the transfer mechanisms and viii chemical species exchange between phases. The finite rate turbulence chemistry interaction model is adopted to demonstrate the strong relation between turbulence and chemistry in reacting turbulent flows. We evaluated the baffles effect on the fluid dynamics in a simplified preneutralizer model, we studied the sensitivity of the reacting flow hydrodynamics to the operating conditions, and we conducted a geometrical optimization study, giving birth to a new preneutralizer design, with optimal gas spargers type and positions. On the other hand, a rotary drum granulator CFD model is tailored. Gas-solid flow is numerically simulated using the Euler-Euler approach, coupled to the Kinetic Theory of Granular Flow (KTGF). The effect of fertilizer particles size on their dynamics is evaluated. We studied the CFD granulation model sensitivity to the drag model subcomponent of the CFD model, and to the key operating conditions. In light of the present study, new operating conditions combinations are proposed for the industrial plant to validate further by an experimental campaign. All the conducted simulations are fully three dimensional with the industrial scale dimensions.

i DEDICATION I would like to dedicate this work to my lovely parents "O & F", my twin sister "Sanae", who taught me to never give up, and to always believe that my future could be whatever I dreamt it to be. Special thanks to my brother in law "Hicham" for his support. I would like to thank them for their everlasting and unconditional love, endless support, their generosity and enormous patience. Without them this work would not be possible! ..... x production.

Ce travail de recherche est composé de deux études différentes; d'une part, la modélisation et la simulation tridimensionnelles du réacteur de preneutralisation, basées sur l'approche Euler-Euler multiphasique, et l'approche à référence multiple (MRF) qui est appliquée pour considérer le mouvement de rotation de l'agitateur. Le modèle CFD gaz-liquide adapté considère les mécanismes de transfert et les échanges d'espèces chimiques entre les phases. Le modèle d'interaction turbulence-chimie à taux fini est adopté pour montrer la relation forte entre la turbulence et la chimie dans les écoulements turbulents réactifs. Nous avons évalué l'effet des chicanes sur la dynamique des fluides dans le preneutralisateur simplifié, et nous avons étudié la sensibilité de l'hydrodynamique de l'écoulement réactif aux conditions de fonctionnement.Nous avons mené une étude d'optimisation géométrique, conduisant à un nouveau design de preneutralisateur, avec un type et des positions optimales des injecteurs de gaz. D'autre part, un modèle CFD de granulateur à tambour rotatif est développé. L'écoulement gaz-solide est simulé numériquement en utilisant l'approche Euler-Euler, couplée à la théorie cinétique des écoulements granulaires (KTGF). L'effet de la taille des particules d'engrais sur leur dynamique est évalué. Nous avons étudié la sensibilité de la granulation au modèle de traînée, dans le developement d'un modèle CFD, et aux conditions de fonctionnement clès de cette opération unitaire. A la lumière de la présente étude, de nouvelles combinaisons des conditions opératoires sont proposées pour l'installation industrielle afin de les valider par une campagne expérimentale.

Toutes les simulations effectuées sont réalisées en 3D, en respectant les dimensions à l'échelle industrielle.

Mots clès : Engrais Di-Ammonium Phosphate, preneutralisation, granulation, flux de réaction multiphasique, modélisation, simulation et optimisation. 3.9 Evolution at increasing gas flow rate F g of gas cavity structure behind the blades of a Rushton turbine. Original imagine from [START_REF] Aw Nienow | Hydrodynamics of stirred bioreactors[END_REF]. . . . . . . . . . . . .

3.10

Interaction between turbulent eddies and fluid particles: a) interaction with several eddies, b) interaction with outer edge of eddy, c) transport with large eddy d) collision with paired eddies [START_REF] Andersson | Modeling the breakup of fluid particles in turbulent flows[END_REF]

3.15

Pumping number for a baffled and unbaffled chemical reactor agitated by a RT [START_REF] Kresta | Advances in industrial mixing: a companion to the handbook of industrial mixing[END_REF]

8.5

Fertilizer particles velocity distributions predicted using the three drag models, for a rotational speed of 8 rpm and 30% filling degree. . . . . . . . . . . . . . . . .

8.6

Front view of a middle transversal plane of particles axial distribution with different drag models at different simulation time with 30% fill rate α s (particles volume fraction), a domain rotational speed of 8 rpm, and an inclination of 3.57

GENERAL INTRODUCTION

The aim of this chapter is to provide the industrial context of this research work, in order to explain the importance and the direct impact of this study on industrial level. A presentation of the targeted objectives and details of the main scientific challenges faced to accomplish the mission are given. Finally, the manuscript structure will be detailed and fellowed by methodological approach. office and the Foreign Exchange Office highlighted the increase in sales of natural and chemical fertilizers [START_REF] Samreen | Phosphorus Fertilizer: The original and commercial sources. Phosphorus-recovery and Recycling[END_REF]. According to the Progressive farming, phosphorusbased fertilizer was 8% more expensive in April compared to March 2022 [START_REF] Mengoub | The Russia-Ukraine War and Food Security in Morocco[END_REF], as prices continue to rise due to the Ukraine crisis, impacting Russian fertilizer exports, which presented a good opportunity for the Moroccan phosphate extraction and transformation leader "OCP Group", the world's largest exporter of phosphate rock and derivatives, managing around 75% of the world's reserves of Moroccan phosphate rock (Analytica, 2022).

In addition to phosphates, OCP exports in the agricultural and agri-food sector stood at 24.37 billion MAD at the end of March against 21.20 billion MAD during the same period last year (Analytica, 2022). The number represents a 29% increase that covered all sectors, including phosphates and derivatives, the agriculture and agri-food sector, and the textile and leather sector in Morocco [START_REF] Amachraa | Morocco emergence in global value chains: Four exemplary industries[END_REF]. In 2022, Phosphate exports registered the fastest annual growth rate across the board, going from 7.2 billion MAD in 2021 to 14.2 billion MAD a year later. Morocco's overall exports reached 62.4 billion MAD at the end of February 2022, a 27.7% rise from last year's 48.9 billion MAD. At the same time, much of the agricultural soil in African countries is of poor quality, partly as a result of geography, environmental degradation and monoculture farming, which has further depleted its health. Rather than just exporting the raw material, Morocco decided decades ago to become one of the world's largest fertilizer producers. In 2020, its state-owned phosphate miner and fertilizer producer OCP Group had a 54% market share of fertilizer exports to Africa. By raising its ambitions, the country has become a "gatekeeper of the world's food supply", according to the policy center report of the 2022 first trimester [START_REF] Mengoub | The Russia-Ukraine War and Food Security in Morocco[END_REF].

The Fertilizer industry in Morocco plays a vital role in ensuring a sufficient market covering by understanding the needs of soils and crops, with a minimal impact on the environment.

Nevertheless, the large demand of phosphates leads to the depletion of good quality reserves, which presents many challenges to deal with, in order to ensure the same quality of fertilizers as requested by users.

The phosphates quality is classified according to their tricalcium phosphate amount, known as Bone Phosphate of Lime (BPL). Thus, a high-grade phosphate contains 70% or more of BPL, an average grade is between 68% and 70%, a low grade is between 61% and 68% and below 61%, phosphate is considered to be of very low BPL with considerable impurities content [START_REF] Amachraa | Morocco emergence in global value chains: Four exemplary industries[END_REF].

To be proactive and to ensure a response to all soil and plants demands, a multi variation of phosphate fertilizer types have been developed (Triple SuperPhosphate (TSP), Single Super-Phosphate (SSP), and Ammoniated phosphate fertilizers: Mono Ammonium Phosphate (MAP), Di Ammonium Phosphate (DAP), and N-P-K fertilizers category, that presents in it self 42 types, containing a fraction of nitrogen, phosphate and potassium, and some trace elements that can be added to provide supplementary nutrients. OCP Group investigates huge efforts to cover all soil needs on fertilizers and reinforce its brand "Feeding the soil to feed the planet".

Ammonium phosphate fertilizers (MAP, DAP, NPK) represent the essential valuable products at OCP Group, especially the DAP fertilizer; this is the reason behind the last interesting investissements (2021-2022) for expanding the storage warehouses DAP (oui). This project aims to increase the storage capacity of this highly demanded fertilizer, in the industrial complex of Jorf Lasfar by expanding 4 warehouses with a capacity of 50 000 tons each, with an allocation of 563 million MAD investment.

Figure 1.1: Sample of Di Ammonium Phosphate fertilizer affected by caking [START_REF] Becker | Phosphates and phosphoric acid: raw materials, technology, and economics of the wet process[END_REF] Among the interesting pillars of a good fertilizer storage, is to have a final product with specific physico-chemical characteristics (Elmisaoui et al., 2021b;[START_REF] Elkarii | A CFD-based surrogate model for predicting slurry pipe flow pressure drops[END_REF] , to ensure the desired quality of the product, and avoid Caking Phenomenom (cf. Fig. 1.1), which consists of crystallization of water-soluble salts and the formation of bridges between the fertilizer grain surfaces during storage. The surfaces also suffer plastic deformation under pressure, and the reduction of the water vapor pressure in the joint between the new contact surfaces causes the particles to adhere to one another [START_REF] Elmisaoui | Etude de l'Impact de la Composition de la Matrice Phosphate sur le Rendement de la Granulation et le Taux de Murissement de l'Engrais Triple Super Phosphate[END_REF].

Two principal explanations, and possible ways to solve this issue, have been given to this phenomenon:

• The chemical compositions of the raw material, a trail that is being dug by applied chemistry researchers;

• The control of the reacting unit operations (Preneutralization and granulation). It is in this mission that we have been entrusted as a Ph.D. doctoral researcher.

Study objectives

One of the main reasons for using MultiPhysics simulations in chemical process engineering is to gain scientific understanding about the physical behavior of complex flows, and their mathematical descriptions in many industrial cases, with respect to the complex details of the geometry, and design. The objectives of this project can be summarized as follows:

1. To develop a new efficient CFD numerical model of the multiphase reacting preneutralizer (PRN), as well as modeling the physical phenomena occurring when considering the interphase exchange mechanisms; considering the turbulence and species transport in the 4. To assess the granulation sensitivity to the main operating conditions: rotational speed, fill rate, and inclination degree; 5. To perform a parametric and structural optimization study of these DAP granulators and to propose innovative solutions for optimal driving and mastery of final product quality.

6. To contribute to the global modeling and a rigorous simulation of the whole granulation loop integrated in the manufacturing processes of DAP [START_REF] Nielsson | Manual of fertilizer processing[END_REF]. Further more the final CFD models will be employed to a digital twin of the main industrial process;

Scientific challenges

There are several scientific challenges to face in order to accomplish the purpose of this research thesis. Some of the major challenges are associated with the complex design of preneutralizer (geometry, gas injectors, ..). Also, to achieve the model convergence with the constraints of considering the three dimensional domain, and real industrial plant dimensions for both studied geometries (PN & RDG). Another difficulty is associated with the challenge of hardware advanced requirements in terms of storage capacity and huge computational time for each calculation and simulation case.

On other hand, the biggest challenge was to develop a relevant CFD model of a reacting multiphase flow considering the turbulence, and to determine the appropriate Turbulence Chemistry

Interactions model allowing to well represent the hydrodynamics in the PN. Noting that this type of research CFD models is mature in the aerodynamics field, specifically; for the CFD modeling of combustion.

Also, coupling the Eulerian Eulerian CFD model to the Kinetic Theory of Granular Flow model for the granular flow hydrodynamics modeling, and choosing the convenient drag model was a challenging step due to the flow complexity in the RDG, and the high sensitivity of CFD simulations results to these sub-models.

This research study was a great opportunity to develop HPC calculation skills to deal with these challenges.

Methodology

A rigorous methodological approach is adopted in this thesis research work to model and analyze the industrial applications of CFD in the preneutralization and granulation process. The approach is based on the following main steps , as in Figure 1.2.

• A literature review is performed to understand the research area and to gather information on related works published by other researchers before proceeding. The process is then analyzed in detail and the scope and the expected outcome from the numerical study are defined based on the formulated hypothesis. To reach the defined goal, sufficient information is gathered about the process to be able to build the numerical domain and to select the most realistic boundary conditions, and physico-chemical parameters.

• The numerical domain is developed based on the information collected in the previous step. The domain is then discretized using a mesh depending on the chosen CFD approach. Suitable mathematical models are then employed to solve a well-defined problem.

• After obtaining the first solution, post-processing is performed. A mesh grid-sensitivity analysis is carried out to ensure the solution grid-independence in the case of mesh-based solvers. The final solution is then validated using available theories, experiments or measurements. If the results are not valid, the model is revisited to make possible adjustments and step 2 is repeated.

• The validated results are then analyzed and hypotheses are made based on the simulation results. Detailed insight regarding the process are provided and possible improvements are suggested in the basis of the phenomena revealed by the CFD simulations.

Manuscript structure

This dissertation is composed of 3 parts, 8 chapters and 2 appendices. It is organized as follows:

• After this general introduction, the first part presents the state of the art of all the fields that have been covered during this work. It is subdivided into three chapters:

1. General Introduction -Chapter 2: DAP chemical fertilizer manufacturing process is presented; every unit operation is described and a special attention is given to the preneutralization and granulation.

-Chapter 3: A literature overview concerning the theoretical aspects of the sparged stirred reactors, and CFD modeling of the hydrodynamics, mass transfer of reacting flows. Some details on complex flows modeling approaches, turbulence modeling and Turbulence Chemistry Interaction (TCI) models are given.

- In light of the present work, new operating conditions combinations are proposed for the industrial plant. An experimental campaign is planned on the industrial plant to further evaluate and compare the new candidates operating points.

• At last, the conclusions of the thesis including the outcomes are presented in chapter 9.

Introduction

A wide variety of phosphate fertilizers exist, allowing to adapt the fertilizer composition to the quality of the soil according to the nutrient requirements. Ammonium fertilizers occupy a large section in the phosphate high value products manufacturing. The chemical processes differ between ammonium (DAP, MAP, NPK, etc...) and non-ammonium (SST, TSP, etc...) based products. In this chapter, we will describe the different components involved in the manufacturing process of Di-Ammonium Phosphate fertilizer. The aim is to describe, present, and investigate the operation units in the chemical transformation loop, with all the ensured functionalities. A special focus will be given to the preneultralization reactor and the rotating drum granulator that constitute the subject of this research work.

Literature review

The world's population increases exponentially and it is expected to reach 10.9 billion by the end of the century. Consequently, food security should be warranted by all the available tools.

Industrials, farmers, and decision-makers are all engaged to enhance food production with application of more environmentally friendly approaches [START_REF] Amit | Fertilizers and Food Production[END_REF][START_REF] Nielsson | Manual of fertilizer processing[END_REF]. To increase crop yields, the demand on chemical fertilizers is a principal pillar in terms of solutions, especially fertilizers compounds based on the three essential elements that are: Nitrogen, Phosphorus, and Potassium (Paul and [START_REF] Victor | Handbook of industrial mixing: Science and Practice[END_REF]. DAP is the leading type of ammonium fertilizer due to its easy applicability in irrigation systems and in the production of suspension fertilizers, since it is totally soluble in water. Furthermore, the fire suppressing applicability of DAP fertilizer in extinguishing forest fire provide necessary nutrients for enhancing trees growth rather than other fertilizers [START_REF] Nielsson | Manual of fertilizer processing[END_REF].

Each year, the fertilizer industry is able to convert millions of tons of phosphate raw materials into a large variety of fertilizers that provide plants the necessary elements for its good growth conditions. These various applications require chemical processes with high precision in terms of control and monitoring of operating conditions (temperature, pressure, stirring contents of chemical elements, level and flow homogeneity).

In 1835, the original Dorrco slurry process was used with a blunger in which a fluid with phosphate salts are mixed in a flash mixer for accomplishing production of the slurry [START_REF] Ivell | Phosphate Fertilizer Production-From the 1830's to 2011 and Beyond[END_REF]. This original process is based for the first time on a particle growth granulation mechanism occurring as a result of application of successive layers from the produced slurry to recycle particles until they achieve the required size. Several improvements have been made to face the different encountered problems with production of patents that have been issued. In 1954, the original Dorr-Oliver process for DAP production began. With three serial ammoniation vessels in which the complete neutralization was carried out. Each stirred reactor operates in a specific Nitrogen-Phosphate rate, with 0.6, 1.4, and 1.85 mole Ammoniation ratio, respectively [START_REF] Nielsson | Manual of fertilizer processing[END_REF].

A special break-through came in 1956 with the Tennessee Valley Authority (TVA) company when Frank Nielsson patented the rotating drum Ammoniator-Granulator instead of the blunger in a new TVA process. This invention showed considerable efficiency and a great impact on the process yield improvement. In the same TVA process, the stirred tank "preneutralizer" vessel was used for partial ammoniation of the phosphoric acid and completion of ammoniation. These modern plants are still the most adopted in the ammonium fertilizers industry, with historical modifications on the designs, operating conditions, and throughput scheming until now [START_REF] Wesenberg | Diammonium phosphate plants and processes[END_REF][START_REF] Kossir | Method for manufacturing an ammonium phosphate fertilizer having low cadmium content, may[END_REF]. The flowchart of the DAP manufacturing process adopted in OCP Group plants is depicted in Fig. 2.1. This process have many benefits: 5. Energy savings, so a reduced industrial OPEX [START_REF] Leikam | Phosphate fertilizers: Production, characteristics, and technologies[END_REF].

As mentioned in Fig. 2.1, the chemical process used at OCP Group for the production of phosphoric acid (ACP) is based on the attack of phosphate rock by sulfuric acid, this type of process is called "wet process manufacturing" [START_REF] Leyshon | Innovations in Slurry Process Granulation Plants[END_REF]Elmisaoui et al., 2021b). At this stage, the granulometry of the phosphate ore plays a crucial role in defining the quality of the attack. This is the reason why phosphate crushing is a conditioning operation that controls the process. Indeed, a high fineness is highly energy consuming, as the latter generates a loss of P 2 O 5 by recrystallization. And this would justify the need for a compromise allowing to obtain an optimal solution.

Manufacturing process description

The flow process diagram of the DAP manufacturing process installed at Jorf Lasfar chemical plants in presented in Fig. 2.2. The transformation loop is composed by five principal operation units from the preneutralizer feeding to the final product [START_REF] Dennis | The Ultimate Granulation Plant[END_REF]. The pipe reactor at the rotating drum is not considered as a functional reactor as assumed for an alternative operation of the preneutralizer stirred tank reactor. This mode of use is not treated in our study, and then all raw material inputs to the pipe reactor are not considered (cf.

Raw material supply

The crude phosphoric acid (54% P 2 O 5 ) is distributed to the various production lines from storage tanks. They are then sent to the various consumption points on each line, namely Phosphoric acid 54% P 2 O 5 and 29% P 2 O 5 are feeding the preneutralizer. Liquid ammonia is pumped from the port to the two huge Ammonia spherical tanks where it is stored in liquid form at -5 • C. Each production line has its own pump to feed the steam heaters, cooler, and granulator.

Sulfuric Acid (ACS) can be considered also as a raw material to produce the DAP. It is used to calibrate the pH level in the preneutralizer.

Preneutralization reaction

The first step of the DAP production is the reaction. For DAP and MAP phosphate fertilizer products, the preneutralization reaction occurs in the preneutralizer (PN) stirred chemical reactor. The PN is a dual diameter stirred tank reactor with a non standard geometry as presented in Fig. 2.3. The smaller diameter is of the bottom part, this specific design ensures a low residence time which enhances species transport and inhibits secondary reactions that fight against P 2 O 5 availability [START_REF] Dennis | The Ultimate Granulation Plant[END_REF]. The upper zone has a large diameter to promote produced gas disengagement. The overall details of the geometrical characteristics of the preneutralizer chemical reactor are given in Fig. 2.3.

Figure 2.3: Preneutralizer chemical reactor

Phosphoric Acid feeds the preneutralizer from the top, and the Ammonia gas is sparged in the bottom via gas injectors with uniform diameter. The PN is equipped by a one agitation stage, that is Pitch Blade Turbine (PBT) located near the bottom, which enhances Ammonia gas dispersion in the useful area. This agitation allows with respect to,other options the following:

• Reduction of foam;

• Good chemical species transport;

• Improvement of the ammonia absorption;

• Homogenization of the reacting flow and temperature within the vessel;

• The suspension of slurry crystals when produced in the PN.

As inter-connectors, two parallel pumps, equipped with independent piping systems, provide the ammonium phosphate slurry to feed the rotary drum from the PN outlet. This arrangement allows the washing and maintenance of one complete production process while the other is operational. This important feature ensures a high duty cycle. The slurry pumps are controlled by valve pumps.

Granulation

The Granulator is the second reacting item in the DAP process and is arguably the most critical unit as it ensures the granules form of the product. Therefore, it has a direct impact on the production quality and capacity [START_REF] Balliu | Performance assessment and model validation for an industrial granulation circuit[END_REF]. Phosphate fertilizer granulation is ensured by rotary drum Ammoniator-granulator that transforms the preneutralized slurry into a granular product of required composition and size [START_REF] Wesenberg | Diammonium phosphate plants and processes[END_REF]. When the N/P mole ratio is equal to 1.0 the MAP is produced, and the preneutralized slurry (MAP) fed to the rotary drum granulator (RDG), on the recycled dry granular material, via the injection pipe [START_REF] Hignett | Fertilizer manual[END_REF].

Figure 2.4: Rotary drum Ammoniator-granulator

The rotational motion of the RDG is assumed to ensure the distribution of slurry liquid on recycled granules surface. Ammonia is sparged by four nozzles through the bottom of AGR, then MAP slurry reacts with Ammonia to form DAP. When the N/P ratio achieved 2.0, the MAP is totally converted to DAP, and the liquid phase is well distributed on layered granules which ensure their size enlargement. A snapshot of the industrial RDG of the DAP fertilizer plant in Jorf Lasfar (107A production line) is shown in Fig. 2.4. It consists of a rotating drum (cylinder)

with a speed of 8 to 10 revolutions per minute (rpm), on six rollers, with four rotational and two concentric, and a toothed crown that serves for rotation. It is also equipped with a kinematic chain composed of a main motor, a reducer and a turner [START_REF] Hignett | Fertilizer manual[END_REF]. The RDG is inclined at 3.57 • towards the outlet in order to facilitate the transfer of a large circulating granular material.

The depth of the granular bed can be adjusted by means of an adjustment at the outlet. The product is conveyed to the dryer through a specially designed chute to minimize clogging and avoid possible damage caused by contact with the hot gas pipes coming from the combustion furnace.

Drying

The dryer allows fertilizer granules to be dried in a co-current air flow. It is dedicated to heat, mass transfer mechanisms and solid transportation efficiency. Drying equipment is composed of two main parts: The furnace fed by fuel liquid and air flow to heat, and the drying zone; that is a long cylinder equipped with shells, where the granular flow meets the hot air flow as presented in Fig. 2.5. Drying zone is equipped with product feed spires at the inlet and lifting shovels to convey the product to the outlet and ensure excellent contact between the dried granules and the hot gases from the combustion furnace [START_REF] Silva | Modelling of fertilizer drying in a rotary dryer: parametric sensitivity analysis[END_REF].

Figure 2.5: Dryer

It is also equipped with five impact hammers to avoid clogging of the first dryer zone. The inside of the shell is equipped with many flights to enhance granules cascading across the full diameter of the dryer and solid contact through the hot gas steam (Arruda et al., 2009a).

The main drying crucial parameter is the residence time required to evaporate the maximum average of moisture from the product, which essentially depends on the granulation mechanisms (agglomerations (15 -30 minutes), layering (5 -10 minutes), ..., etc.) that allows enough residence time to diffuse moisture inside the solid, through pores, to the surface. Generally, the MAP and DAP fertilizer granules drying process takes 5 to 7 minutes [START_REF] Kossir | Method for manufacturing an ammonium phosphate fertilizer having low cadmium content, may[END_REF].

As a rotary drum, the dryer participates indirectly to the increase of process yield, with a good control of its rotational speed, its slope (around 1.5 • and the gas flow temperature). For Ammonium phosphate fertilizers, the maximum solids temperature exiting the dryer is about 90 • C (significantly less for NPK's containing urea) [START_REF] Dennis | The Ultimate Granulation Plant[END_REF]. Any exhaust temperature can generate ammonia fumes more than the high amount of energy dissipation [START_REF] Nielsson | Manual of fertilizer processing[END_REF].

Classification: Screening & Crushing

The classification stage in the manufacturing process is composed by two successive unit operations, that are: screening and crushing.

A key element in the manufacturing of DAP fertilizer production is screening, which ensures the high precision in granules size and physical characteristics of the product. However, a great screening system by well designed screens allows to avoid a supplementary downstream screening of the final product (polishing before storage) [START_REF] Keßler | Modern Screening Technologies for Maximizing Quality & Efficiency Improvements[END_REF]. DAP screening unit is composed of four double deck screens equipped with vibrating motors. All the recycling products pass through these screens.

The dryer elevator feeds a screen feed conveyor which transports the recycled material to a system of chutes and flap boxes for distribution on the four screens [START_REF] Ivell | The evolution of screening systems for optimum granular fertilizer product quality[END_REF].

As depicted in Fig. 2.6, to improve screening efficiency, each screen is equipped with a vibrating feeder to ensure product distribution across the entire screen width. The screened fine fertilizer granules are conveyed to the recycling belt by means of chutes [START_REF] Keßler | Modern Screening Technologies for Maximizing Quality & Efficiency Improvements[END_REF], the coarse granules rejected from upper screen, is conveyed to the crushing unit, and on-size DAP granules with diameter from 2 mm to 5 mm are moved to the conditioning unit. The final product ratio is adjusted by the shutter box and the amount of recycled product in the classification system, which is an automated process [START_REF] Schultz | Production of granular NPKs in ammonium phosphate plants: some important differences[END_REF].

Furthermore, the quality for these final products is often expected to be more than 90% meaning that only 10% of the granules may be outside of the specified range. The complementary facility in the classification stage is crushing, which ensures size reduction of the oversize fertilizer particles. It is composed by a set of crushers to avoid the risk of recycling a large amount of uncrushed or partially crushed oversized granules to the rotating drum granulator.

This scenario can upset the granulation that is very sensitive to the particle size distribution

Conditioning: Cooling & Coating

The conditioning unit in 107 A manufacturing line of DAP fertilizer is composed by a cooler, a coater, and an intermediate secondary screens as highlighted in Fig. 2.7. The cooler is an industrial fluidized bed unit designed to reduce the product's temperature to approximately 50 • C with ambient air that has been cooled to 27 • C by two air conditioners, which prevents caking, creates a more stable granules, and mitigates the potential for any heat damage to downstream unit operation [START_REF] Nielsson | Manual of fertilizer processing[END_REF]. Each air conditioner is an oil heat exchanger that ensures the ambient air temperature to 21 • C by evaporating a portion of the liquid ammonia fed to the unit.

The saturated air is then heated by steam circulating through a 27 • C coil to reduce its humidity. The two cooler fans draw the air through the two conditioners and discharge it at a slightly higher pressure to the cooler [START_REF] Dennis | The Ultimate Granulation Plant[END_REF]. The product exiting the cooler is discharged to an elevator that feeds the secondary screens for the particles size validation. The secondary screens output feeds the rotary drum coater, in which the coating material (generally an amino oil) is pulverized on fertilizer granules to promote a good physical condition (flowability) during storage and handling. It is a crucial step for some phosphate fertilizers with a high caking tendency. The product at the exit of the coater is discharged on the conveyor of the The conditioning unit (cooler and coater) is a secondary component in fertilizer manufacturing. Its application depends on the subsequent processing or handling required for specific physio-chemical properties of the final product or customer's requirements (oui; Wesenberg, 2018).

Focus of this study

This research work aims to use CFD approach to investigate the hydrodynamic flows in the reacting unit operations of the DAP Manufacturing process. From the industrial context of that research project, the preneutralization and the granulation are the mastering operations, and the challenging phenomenology in terms of adaptability to the raw material quality. Despite several attempts and successful research, understanding the hydrodynamics in the preneutralizer and the rotary drum granulator still cause major challenges. The complex hydrodynamics and fluid dynamics involved in these unit operations make their monitoring really challenging. Furthermore, the sensitivity of both operations to the operating conditions and geometrical characteristics are studied with special attention in the present study by gradually increasing the complexity of systems. The objective of this thesis is to be able to understand the hydrodynamics of multiphase flows in the PN and RDG from the results of that work, which is a prerequisite for the development and the intensification of these unit operations to be easily adaptable to the raw material quality changes.

Conclusions

The fundamentals of particle formation and processing for agglomerated Ammonium Phosphate fertilizers differ markedly from other fertilizer categories. For this reason, many welldesigned and understood manufacturing plants experience difficulties of monitoring. The difficulties are most pronounced when raw material quality changes considerably. The plant design and process components described in this chapter are the principal unit operations involved in the production of DAP at OCP Group plant at Jorf Lasfar. This chapter highlighted also the reason behind the choice of the preneutralization and the granulation to study in this thesis research project.

In the next chapter, we will present the general theoretical aspects regarding the multiphase systems: stirred tanks and rotary drums. The main notions about hydrodynamics, the mass transfer occurring this units and also their numerical study through CFD approach will be presented, to provide the key concepts required to understand the contents presented throughout the manuscript.

CHAPTER 3

Introduction

The overall objective of that chapter is to gain insight into the hydrodynamics prevailing in the sparged stirred preneutralizer reactor. In this chapter, the theoretical background employed to analyze the principal aspects of the turbulent industrial mixing systems hydrodynamics, industrial gas-liquid stirred reactors are first introduced. Afterwards, the governing equations of fluids mechanics and turbulence models, especially the Reynolds Average Navier-Stokes Simulation (RANS) turbulence modeling are described, and the different modeling assumptions required to handle turbulent reacting flows are discussed. The existing mathematical models and numerical schemes used to model this complex domains are given, with details on inter-phase mass transfer mechanisms, transport phenomena modeling and the mixing systems hydrodynamics evaluation parameters. Finally, the applicability of this knowledge to the process of interest is investigated.

Theoretical aspects on industrial mixing

A variety of techniques, including agitation, sparging, and static flow manipulation, can be used to enhance mixing in industrial mixing processes. One of the most frequent mixing unit operations is the stirred tank reactor (STR), used for the first time by Georg Bauer in 1556 for mining and metals treatments [START_REF] Nienow | Stirring and stirred-tank reactors[END_REF]. Despite, the large use of stirred tank reactors at industrial scales instead of other reactors types (bubble column, packed bed, plug flow), they still attract significant attention in terms of engineering research, and development, with interesting paths of innovative intensification and optimization. Also STRs present one of the biggest challenges for computer simulation [START_REF] Kresta | Advances in industrial mixing: a companion to the handbook of industrial mixing[END_REF]. Stirred tanks normally have an impeller installed on a shaft, with additional interior components such as coils, or baffles as an alternative. Computational fluid dynamics (CFD) must take into account a variety of factors while modeling a stirred tank (Paul and [START_REF] Victor | Handbook of industrial mixing: Science and Practice[END_REF].

From many years, the understanding of mixing in Stirred Tank Reactors (STRs) progressed significantly [START_REF] Mann | Computational fluid mixing for stirred vessels: progress from seeing to believing[END_REF]. The first attempt to bring engineering into the mixing art was done by James Thompson in 1855 [START_REF] Nienow | Stirring and stirred-tank reactors[END_REF]. He measured the power required to rotate a disc in water and found a relation that is still used to predict the power consumption of STRs (Paul et al., 2004). Stirred tank reactors (STRs) are used with different liquids categories, slurries, very low gas flow rates, and large liquid volumes [START_REF] Charpentier | Mass-transfer rates in gas-liquid absorbers and reactors[END_REF], also are very known by the well-mixed state achievement easily, which aids in providing necessary hydrodynamics requirements [START_REF] Mcfarlane | Studies of high solidity ratio hydrofoil impellers for aerated bioreactors. 2. Air-water studies[END_REF]Elmisaoui et al., 2021b).

Typical and standard STR units have a small height-to-diameter ratio relative to other reactor types, where the diameter D can vary from about 0.1 m for experimental units to 10 m for industrial applications. The impeller and baffle dimensions, as well as the impeller clearance are typically a specified fraction of the tank diameter [START_REF] Charpentier | Mass-transfer rates in gas-liquid absorbers and reactors[END_REF].

In some specific applications of STRs this rules of sizing are not respected, especially for industrial multiple impeller designs that are mostly dependent to many practical considerations [START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF].

CFD approach for complex flows modeling

Computational Fluid Dynamics (CFD) is first used in industrial manufacturing to complement experimental data for obtaining results either experimentally difficult to deduce from sensitivity studies that would be time-consuming experimentally [START_REF] Essam | CFD history and applications[END_REF]. The objective behind is to develop a better comprehension of physical phenomena and to provide further understanding of what happens in chemical reactors. It is also used for designing equipment which has to be linked with Computer-Aided Design (CAD) use in the mechanical industry. More recently, it has been considered as an efficient tool to be used before performing any experiments and of great help for designing the corresponding experimental set-ups.

The first step in CFD modeling is called "pre-processing". It consists of setting up the case, defining the geometry and most importantly, since it is less visible, choosing the adequate boundary conditions and the appropriate mesh. With the increasing Central Processing Unit (CPU) performances and some novel meshing approaches, young researchers tend to pay less attention to the mesh quality [START_REF] Kramer | Modeling of industrial crystallizers for control and design purposes[END_REF]. However, it is well known that a carefully prepared mesh with low skewness factor cells and oriented in the main flow direction will contribute to a more efficient convergence, and that minimum numbers of cells are required in zones with high shear. Boundary condition choice is of even greater importance since it gives the physical meanings to the calculations. The choice is often not obvious since the flow characteristics at inlets (phase distributions, particle sizes, pressure, etc.) are not always well known. Thus can be the subject of sensitivity analysis during the calculation step. Thus, this pre-processing step is the period of time in which researchers invest to make the real development phase more efficient.

The second step "Processing"; consists of the calculations themselves, which corresponds, in the present analogy, to the process development step. This is the hard part, concentrating most of the work. In CFD simulations, appropriate models (in terms of physical modeling and mathematical solvers) will have to be chosen. Calculations will be performed directly or via a multiscale approach so as to develop better understanding and physical correlations that could be used in process flow sheet simulations or in process design rules [START_REF] Bartels | Computational fluid dynamics applications on parallel-vector computers: computations of stirred vessel flows[END_REF]. The purpose of this step is to develop the appropriate knowledge to design the chemical reactors, and especially the reactive section, and further bring corresponding innovations to manufacturing process. The third step is called "post-processing" or post-treatment. It consists of making the most of the calculations; results data treatment and either to show colourful illustrations or provide quantitative information that would increase the value of simulations. Concerning the process development, the analogy would correspond to diagnostics, troubleshooting or equipment design until achieving the desired yields or performances.

CFD is able to predict fluid flows, chemical reaction rates, mass and heat transfer rates, and other occurrences by solving a set of appropriate mathematical equations [START_REF] César | Three-dimensional CFD simulation of a large-scale grate-fired biomass furnace[END_REF][START_REF] Md | Literature review of accelerated CFD simulation methods towards online application[END_REF]. CFD also provides useful information for regions with intense or mild turbulence zones, Reynolds stresses, vortex structures, circulation patterns, flow behavior, and many other parameters. Not with standing its potential, fundamental knowledge in multiphase flow simulations is still lacking because of its complexities, particularly in stirred vessels that have moving components [START_REF] Echi | CFD simulation and optimization of industrial boiler[END_REF][START_REF] Saad Benjelloun | Multi-fluid flows: mathematical analysis, modelling and simulation[END_REF].

We can summarize the CFD interests for multiphase flows modeling in the following points:

• CFD can be used to augment design correlations and experimental data.

• CFD provides comprehensive data that are not easily obtainable from experimental tests.

• CFD reduces scale-up problems, because the models are based on fundamental physics and are scale independent.

• When evaluating industrial scale problems, CFD is often used to understand the root cause of problems, not just the effects.

• CFD can be used to complement physical modeling. Many research works actually uses CFD to analyze new systems before deciding which and how many validation tests need to be performed. Finite Volume Method-based CFD solvers are particularly able to deal with complex multiphase flow and heat transfer problems. They are mesh-based solvers that have been used extensively for many years to resolve immensely difficult problems [START_REF] Lim | Computational Fluid Dynamics[END_REF]. Nevertheless, a significant difficulty for such solvers and the general method continues to be the solution of industrial issues requiring high deformation, numerous solid motions, strong and complicated fluid-solid or fluid-fluid interactions. Industrial mesh creation difficulties involving fluid-fluid or fluid-solid interactions need for special caution [START_REF] Gershtein | Computational fluid dynamics in industrial combustion[END_REF]. In the literature, the limitations, difficulties, and challenges in resolving complex industrial problems have been thoroughly examined [START_REF] Frank | Advances in computational fluid dynamics (cfd) of 3-dimensional gas-liquid multiphase flows[END_REF].

Multiphase flows are encountered in many processes in industrial operations. Numerous examples can be found in the chemical, petroleum, pharmaceutical, agricultural, biochemical, food, electronic, and power-generation industries. Due to the inherent complexity of multiphase flows, from a physical as well a numerical point of view, "general" applicable computational fluid dynamics (CFD) codes are not available [START_REF] Van Wachem | Methods for multiphase computational fluid dynamics[END_REF].

Reasons behind the lack of fundamental knowledge on multiphase flows are:

1. Multiphase flow is a very complex physical phenomenon where many flow types can occur (gas-solid, gas-liquid, liquid-liquid, etc.) and within each flow type several possible flow regimes can exist (annular flow, jet flow,bubbly flow, etc.) [START_REF] Rielly | Computational fluid mixing. Food Mixing: Principles and Applications[END_REF].

2. The complex physical laws and mathematical treatment of phenomena occurring in the presence of the two phases (interface dynamics, coalescence, break-up, drag, ... ) are still largely undeveloped. For example, to date there is still no agreement on the governing equations.

3. The numerical methods for solving equations and closure laws of multiphase flows are extremely complex. Very often multiphase flows show inherent oscillatory behavior, requiring costly transient solution algorithms. Almost all CFD codes apply extensions of single-phase solving procedures, leading to diffusive or unstable solutions, and require very short time-steps [START_REF] Marchisio | Multiphase reacting flows: modelling and simulation[END_REF].

Impeller rotational movement modeling

As the impeller is the source of movement in a STR's turbulent flow, its precise representation is crucial to the flow and mixing solution. Several approaches were employed in early studies to provide the impeller's boundary conditions, such as experimental data, momentum sources, etc. Considerable multi zone solution approaches are available [START_REF] Vivek V Ranade | Computational flow modeling for chemical reactor engineering[END_REF]. The famous two main approaches for the simulation of impeller rotation in a mixing vessel are: Sliding mesh approach (SM), and Multiple Reference Frame (MRF) approach [START_REF] Hadane | Cfd investigation of the agitation in the desupersaturation during the wet-process phosphoric acid (wppa) process[END_REF].

The principle of multi zone resolution methods, is that the flow around the impeller is solved using such a reference frame that rotates with the impeller. In this frame, the impeller becomes a fixed object, avoiding the need for a moving grid (Paul et al., 2004). During the solution procedure, the interface boundary conditions between the revolving and stationary frames are matched. When RANS models are employed for turbulence modeling, most reported findings demonstrate that the sliding mesh approach has relatively little benefit over the quasi-steady MRF method for both mean flow prediction and space-averaged hydrodynamic parameters calculations such as power number and flow number [START_REF] Kresta | Advances in industrial mixing: a companion to the handbook of industrial mixing[END_REF].

The Multiple Reference Frame (MRF) method allows modeling of agitated baffled tanks of different complexities. A rotating frame (i.e. co-ordinate system) is used for the zone containing the rotating components while a stationary frame is for stationary regions (cf. Figure 3.1). Inside the rotating frame, momentum equations are solved in the frame of the enclosed impeller, while those outside are solved in the stationary frame. As the solution proceeds, the MRF interface continuously transfers information. [START_REF] Cullen | Food mixing: Principles and applications[END_REF] We suppose the rotating frame at position r 0 relative to the stationary frame as shown in transformed from the stationary frame to the rotating frame using the following formulas:

v r = v -u r , (3.1)
where, v r is the relative velocity seen from the rotating frame, v is the absolute velocity seen from the stationary frame, and u r is the "whirl" velocity due to the moving frame given by u r = ω × r.

(3.2)

Figure 3.2: Fixed and rotational reference frames

When the equation of motion is transferred to the rotating reference frame, the continuity and momentum equations respectively become (3.4) where; (2ω × v r ) and (ω × ω × r) respectively are the Coriolis and centripetal accelerations, and τ r is the stress tensor based on v r . The momentum equation for the absolute velocity is

∂ρ ∂t = -∇ • ρv r , (3.3) and ∂ (ρv r ) ∂t = -∇ • (ρv r v r ) -ρ (2ω × v r + ω × ω × r) -∇p + ∇ • τ r + ρg + F ,
∂(ρv) ∂t = -∇ • (ρv r v) -ρ(ω × v) -∇p + ∇ • τ + ρg + F (3.5)
where ω × v embodies the Coriolis and centripetal accelerations. The MRF method is recommended for industrial scale simulations of stirred vessels.

In this research study, MRF approach is adopted to model the turbine motion in the preneutralizer reactor.

Governing equations and mathematical models

The governing equations upon which every CFD model is based on are the Navier-Stokes equations, which represent the local conservation of mass and momentum, respectively [START_REF] Bird | Transport phenomena[END_REF]:

∂ρ ∂t + ∇ • (ρu) = 0, (3.6) ∂ρu ∂t + ∇ • (ρuu) = -∇P + ∇ • τ + g, (3.7)
where, τ is the viscous stress tensor and g is the vector of the body forces per unit of volume acting on the fluid (gravity force in most cases). The stress tensor is given by

τ = µ ∇u + ∇u T - 2 3 ∇ • (uI) , (3.8)
where I is the identity matrix. It can be seen that the last term in Eq. (3.8) is zero for incompressible fluids. In this case, the stress tensor becomes

τ = µ ∇u + ∇u T = µD, (3.9)
here, D is the rate-of-deformation tensor.

Time-averaging over a reasonable time permits to describe separately fluctuations and nonturbulence (or statistically) quantities. Applying Reynolds decomposition to both pressure and velocity in Eqs. (3.6) and (3.7) gives (3.11) where the strain-stress tensor is expressed for a Newtonian fluid. Time-averaging leads to a new set of PDE known as Reynolds-averaged Navier-Stokes (RANS), presented in Eqs. (3.12) and

∂ρ ∂t + ∇ • (ρu) + ∇ • (ρu ′ ) = 0, (3.10) ∂(ρu) ∂t + ∂ (ρu ′ ) ∂t + ∇ • ρ (u + u ′ ) (u + u ′ ) = -∇ P -∇P ′ + µ ∇ 2 u + ∇ 2 u ′ + g,
(3.13). The time average operation is done over a timescale larger than those of turbulence, but smaller than the timescale of mean quantities. This means that time derivative terms in these equations take into account a variation that takes place during a time larger than turbulence timescales:

∂ρ ∂t + ∇ • (ρu) = 0, (3.12) ∂(ρu) ∂t + ∇ • (ρu) = -∇ P + µ∇ 2 u + g + ∇ • τ ′ , (3.13)
The averaging generates a new term, τ ′ = -ρu ′ u ′ , also referred to as the Reynolds stress tensor (RST). It is an unknown symmetric second-order tensor which has to be modeled in order to close the problem. An approach to do would be to derive a transport equation for the RST itself; unfortunately, this leads to a third-order tensor. Deriving an equation for this term would yield an equation that contains fourth-order moments of the velocity components. This goes on indefinitely and is noted as the closure problem. The Boussinesq approximation states that RST is proportional to the mean velocity gradients, describing turbulent transport of momentum as a diffusivity process: a turbulent viscosity is introduced and its modeling is the aim of the several turbulence models that are based on Boussinesq's hypothesis:

τ ′ = -ρu ′ u ′ = - 2 3 (ρk + µ t ∇ • u) I + µ t ∇u + ∇u T . (3.14)
Here, k is the turbulent kinetic energy, equal to k = 1/2 u ′ • u ′ . Turbulence models are classified as zero-, one-, or two-equation models if they include zero, one, or two PDE to solve in conjunction with the RANS equations Eqs. (3.12) and (3.13). Only two-equation models were employed in this work, hence only these will be discussed further.

Meshing of solution domain

The mesh type has a significant impact on solution accuracy. To discretize a continuous flow domain, several mesh types are available. Refer to Figure 3.3 for a comprehensive list of mesh types. It is generally known that hexahedral (hex) cells produce a more accurate solution than tetrahedral (tet) cells. Tet cells, on the other hand, are frequently considerably easier to design for sophisticated geometry [START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF]. This is especially true in industrial modeling. As a result, it is a trade-off between the simplicity of mesh production and the precision of the solution. Another drawback of tet cells is that the model is sometimes many orders of magnitude bigger than that of hex cells. This increases memory use as well as CPU time. Thus, for a mesh that will be used repeatedly for parametric research or for more accurate answers, the higher initial cost of creating a set of hex meshes frequently pays off in the long simulations run (Paul et al., 2004).

As the grid size decreases, the simulation results become grid-independent. A grid independence analysis is carried out in order to create a grid size-independent solution. Grids of different sized are produced in this study, and simulations are run on each grid size. Grid independence is claimed to have been attained when the solution does not change with grid size.

When the grid size for a particular solution domain is reduced, the number of grid cells increases in order to attain high accuracy at the price of processing time [START_REF] Krishna | Three-phase eulerian simulations of bubble column reactors operating in the churn-turbulent regime: a scale up strategy[END_REF]. In this context, "sufficiently small" means that the mesh width (or grid cell size) must not be too large in order to avoid influencing convergence; on the other hand, too small grid size yields a large number of cells for which the discretized equations must be solved, quickly resulting in a large computational demand in terms of processor time and memory usage: In a threedimensional calculation, reducing grid cell size to half its original value results in an eightfold increase in the number of equations to be solved a number that is proportional to computational demand.

Figure 3.3: Mesh types to use in computational grids development [START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF] After meshing is complete, the mesh's quality must be confirmed. The skewness of each cell must be computed for this reason. Skewness describes how much the cell departs from a conventional optimal cell shape, where the optimum cell size is the same circumradius as an equilateral cell [START_REF] Mardani | Numerical investigation of a double-swirled gas turbine model combustor using a rans approach with different turbulence-chemistry interaction models[END_REF].

The discretization process produces a finite set of coupled algebraic equations that must be solved simultaneously in each cell of the solution domain. Because the equations that govern fluid flow and related processes are nonlinear, an iterative solution approach is necessary.

Numerical resolution procedure

After meshing is complete, simulations must be run. There are several solvers available for this purpose, including FLUENT, CFX, POLYFLOW, and STAR CCM. All solutions fulfill the fundamental requirements of combining fluid characteristics and applying suitable boundary conditions. They do, however, each have distinct advantages and disadvantages. To discretize partial differential equations to ordinary difference equations, the finite volume method is frequently utilized, and numerical methods are used to solve them using the appropriate initial and boundary conditions. and associated to its normal vector. The conservation law for a specific quantity u in the control volume V means that its variation in time ∂ ∂t udV is equal to the contributions of the internal and external sources or applied forces [START_REF] Cullen | Food mixing: Principles and applications[END_REF], .

Figure 3.4: Fixed control volume for conservation laws balances

The control volume integration, which forms the key step of the finite volume method that distinguishes it from all other CFD techniques. For each finite size cell, the resultant statements express the (accurate) conservation of key characteristics. One of the key advantages of the finite volume technique is the obvious relationship between the numerical algorithm and the underlying physical conservation principle, which makes its ideas considerably easier for engineers to understand than the finite element and spectral methods. Within the finite control volume (Figure 3.4), the conservation of a general flow variable, such as a velocity component or enthalpy, may be described as a balance between the many processes that tend to increase or reduce it. We have the following words:

                         Rate of change of φ in the control volume with respect to time                          =                                Net rate of increase of φ due to convection into the control volume                                +                                Net rate of increase of φ due to diffusion into the control volume                                +                          Net rate of creation of φ inside the control volume                          (3.15)
CFD codes include discretization techniques that can be used to treat the key transport phenomena of convection (transport due to fluid flow) and diffusion (transport due to variations in φ from point to point), as well as the source terms (associated with the creation or destruction of φ) and the rate of change with respect to time. The underlying physical phenomena are complicated and non-linear, CFD codes include discretization techniques that can be used to treat the key transport phenomena of convection (transport due to fluid flow) and diffusion (transport due to variations in φ from point to point), as well as the source terms (associated with the creation or destruction of φ) and the rate of change with respect to time.

The velocity field, in general, is unknown and emerges as part of the entire solution process alongside all other flow variables. This section examines the most often used techniques for calculating the whole flow field. The pressure gradient term, which forms the main momentum source term in most flows of engineering importance. If the pressure gradient is known, the procedure for getting discretized equations for velocities from momentum equations is identical to that for any other scalar. So, specific resolution agorithms are used. In general-purpose flow computations we also wish to calculate the pressure field as part of the solution, so its gradient is not normally known beforehand. However, if the flow is incompressible, the density is constant and hence not related to the pressure. In this situation, the pressure-velocity coupling tends to add a constraint in the flow field solution: if the correct pressure field is utilized in the momentum equations, the resulting velocity field should satisfy continuity.

An iterative solution strategy may be used to solve both the difficulties associated with the nonlinearities in the equation set and the pressure-velocity linkage, such as the SIMPLE algorithm. The abbreviation SIMPLE refers for Semi-Implicit Method for Pressure Linked Equations. The algorithm was initially put up by Patankar and Spalding (Versteeg and Malalasekera, 2007) and is simply a guess-and-correct approach for the estimation of pressure on the staggered grid pattern. As described in Figure 3.5. In this algorithm the convective fluxes per unit mass F through cell faces are evaluated from so-called guessed velocity components. Furthermore, a guessed pressure field is used to solve the momentum equations, and a pressure correction equation, deduced from the continuity equation, is solved to obtain a pressure correction field, which is in turn used to update the velocity and pressure fields. To start the iteration process [START_REF] Saeed Hussin | Computational analysis of two-dimensional wing aeroelastic flutter using Navier-Stokes model[END_REF] we use initial guesses for the velocity and pressure fields. As the algorithm proceeds our aim must be progressively to improve these guessed fields. The process is iterated until convergence of the velocity and pressure fields.

Turbulence modeling

Turbulence is a continuous phenomenon that exists on a large range of length and time scales, which are still larger than molecular scales. In order to visualize turbulent flows, one often refers to turbulent eddies, which can be thought of as a local swirling motion whose characteristic dimension is on the order of the local turbulence length scale. Turbulent eddies also overlap in space, where larger eddies carry smaller ones. As there exists a large range of different scales (or turbulent eddy sizes), an energy cascade exists by which energy is transferred from the larger scales to the smaller scales, and eventually to the smallest scales where the energy is dissipated into heat by molecular viscosity. Turbulent flows are thus always dissipative.

Turbulent flows also exhibit a largely enhanced diffusivity. This turbulent diffusion greatly enhances the transfer of mass, momentum, and energy. The apparent stresses, therefore, may be of several orders of magnitude greater than in the corresponding laminar case.

There are several numerical methods for solving the governing equations and dealing with turbulent flows. The Reynolds averaged Navier-Stokes (RANS) method begins with ensemble averaging the governing equations, which employs the Reynolds decomposition. The turbulent stresses are used to explain the turbulence, and information on the averaged quantities is provided. Large eddy simulation (LES) is an intermediate tool for computing turbulent flows. In LES, the large scales of turbulent flow are resolved and the small scales of the flow are modeled, using the information obtained from the resolved field. Unlike RANS techniques, LES provides a more realistic and detailed time-resolved data analysis possibility, which gives some knowledge about different structures of turbulent flows. [START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF] Indeed, LES is a highly effective technique for turbulent flow calculations, especially when no response is involved and modeling difficulties are limited to turbulence at small scales. The most accurate computational method is direct numerical simulation (DNS), which solves all flow equations without making any additional assumptions and resolves all scales from large integral scales to small Kolmogorov scales of the flow and even fine structures of the flame. The time and spatial resolution requirements for the DNS in turbulent combustion flows are more rigorous than in non-reacting counterpart flows.

RANS based turbulent models

In this study we had choose a two equation turbulence models for modeling turbulence, as The strategy adopted in RANS method is to set aside the instantaneous fluctuations of the flow, in order to express the mean field equations. The problem of simulating all the small structures of the flow is thus eliminated. Mesh size is not required to be very small as in DNS and LES. This approach is therefore less expensive in terms of calculation time than the DNS and the LES, which explains its wide use in the industrial scale CFD models. Nevertheless, the RANS approach is empiric, which makes it unreliable in some configurations. It consists on solving the averaged Navier-Stokes equations obtained by introducing a decomposition of the variables of the problem into a mean value and a fluctuation term.

In kε models, the turbulent length scale is defined in terms of turbulent dissipation rate ε. The standard formulation is then based on the two transport equations for k and ε; for an isothermal and incompressible system, these read:

∂(ρk) ∂t + ∇ • (ρku) = ∇ • µ + µ t σ k ∇k + G k -ρε, (3.16) ∂(ρε) ∂t + ∇ • (ρεu) = ∇ • µ + µ t σ ε ∇ε + C 1ε ε k G k -C 2ε ρ ε 2 k , (3.17)
G k is the production of turbulent kinetic energy due to mean velocity gradients, and µ t is the turbulent viscosity; these quantities are defined as:

G k = -ρu ′ u ′ : (∇u)µ t = C µ ρk 2 ε (3.18)
In the model, there are five constants, the values of which are:

C 1ε = 1.44, C 2ε = 1.92, σ k = 1.0, σ ε = 1.3, and C µ = 0.09.
The realizable kε model, is based on a different transport equation for the turbulent dissipation rate. Moreover, C µ is not a constant anymore, but a function of the turbulence quantities and the velocity field. The following equation is solved to obtain ε :

∂(ρε) ∂t + ∇ • (ρεu) = ∇ • µ + µ t σ ε ∇ε + ρC 1 γε -ρC 2 ε 2 k + √ νε , (3.19)
where, C 1 = max 0.43, η(η + 5) -1 , η = γk/ε, and γ = √ 2S : S. S is the rate-of-strain tensor, which represents the symmetric part of the rate-of-deformation tensor D (i.e. S = 1/2D ).

The model constants are: C 2 = 1.9, σ k = 1.0, and σ ε = 1.2.

The realizable kε model has been extensively validated for a wide range of flows and its performance has been found to be substantially better than that of the standard kε model.

Multiphase flow model formulation

Multiphase stirred reactors often operate in a dispersed flow regime, with dispersed bubbles and/or solid particles in the continuous phase. Dispersed and segregated (or stratified) multiphase flows are the two major classifications of multiphase flows. If good results can still be obtained when the dispersed phase is constituted of small particles (such as bubbles or droplets), as either the volume fraction of the dispersed phase or the particle dimension increases, the dispersed phase heavily influences the behavior of the continuous phase, which increases the simulation difficulties [START_REF] Mcfarlane | Studies of high solidity ratio hydrofoil impellers for aerated bioreactors. 2. Air-water studies[END_REF].

For multiphase dispersed systems, there are essentially two models: the Euler-Lagrange model and the Euler-Euler model. The frame in which the flow governing equations are formulated differs between these formulations. The continuous medium is described in an Eulerian frame of reference (the point of view is stationary), whereas the dispersed phases are described in a Lagrangian frame of reference. This approach produces excellent results, allowing one to describe the motion of each particle, but it is plainly unsuited for systems with large levels of gas holdup, due to the huge computing cost involved.

On the other hand, the Euler-Euler approach, interprets the phases as interpenetrating continuous media, solving continuity and momentum equations for each phase. The disadvantage is that closure equations are required to characterize the interaction between the flow phases.

Due to the fraction of gas to liquid phase, Euler-Euler approach will be considered in that work of preneutralization CFD modeling.

The Euler-Euler model, also known as the two-fluid model, is obtained by volume averaging the governing equations, taking into account the bulk properties of each fluid, and introducing the local volume fraction k, where k represents the k-th phase (e.g., k = 1 for liquid phase and k = 2 for bubbles). As a result, Eqs. (3.6) and (3.7) become:

∂α k ρ k ∂t + ∇ • (α k ρ k u k ) = 0, (3.20) ∂α k ρ k u k ∂t + ∇ • (α k ρ k u k u k ) = -∇ (α k P ) + ∇ • (α k T k ) + α k g k + f k , (3.21)
where f k is the interaction force between the k-th phase and the others. Gas-liquid multiphase flow systems in chemical engineering accounts for about 25% of all chemical reactions [START_REF] Beenackers | Mass transfer in gas-liquid slurry reactors[END_REF]. Gas dispersion requires that the entire reactor volume be used to mix the gas into the liquid. The low solubility of most gases limits gas absorption to the point that gas-liquid mass transfer becomes the rate limiting step for the overall reaction [START_REF] Bouaifi | A comparative study of gas hold-up, bubble size, interfacial area and mass transfer coefficients in stirred gas-liquid reactors and bubble columns[END_REF]. This limitation is even more severe in systems using very low solubility gases. Some of which are crucial in industrial applications. Thus, increasing the gas-liquid mass transfer is the fastest strategy to boost productivity in these chemical processes. [START_REF] Gary B Tatterson | Fluid mixing and gas dispersion in agitated tanks[END_REF].

The mass transfer between the gas and the liquid phase is usually described by considering the two-film theory of Lewis-Whitman (Kutateladze and Nakoriakov, 1984). It consists on describing the exchange of a species A from a specific phase to another, because of the presence of a concentration gradient. This driving force is focused at the interface, within each phase's film 3.7.

Given the component A being transferred from the gas phase to the liquid phase, the flux for each phase can be written as law [START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF].

J A,g = k G p A -p A,i , (3.24) J A,l = k L C A,i -C A , (3.25)
p A,i = H A C A,i , (3.26)
where H A is the Henry constant for component A at a specific temperature and for the considered liquid. Introducing the equilibrium constant, Eq. (2.24) becomes

J A,g = H A k G C * A -C A,i , (3.27)
here C * A is the equilibrium concentration in the liquid. Dealing with concentrations at the interface is not easy, so an overall mass transfer is usually introduced, so the flux can be expressed in terms of the overall driving force respect to the liquid concentration:

J A,l = K L C * A -C A .
(3.28)

Rearranging the previous equations, and considering the equality of the fluxes through each phase at the steady state, the different mass transfer coefficients can be related as follows

1 K L = 1 k L + 1 H A k G ≈ 1 k L . (3.29)
There are several methods for expressing k L , but the most frequent are based on Higbie's penetration theory and the surface renewable model [START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF].

Empirical correlation is commonly used to evaluate the volumetric mass transfer coefficient in stirred tanks. These correlations represent k L a as a function of operating circumstances and take the following general form

k L a = P g V a U b g µ c .
(3.30)

By fitting experimental findings, many writers offered alternative values for the exponents a, b, c. This large range of values can be explained by the fact that all of the authors used various measuring techniques and conducted tests at varied gas velocities, tank diameters, and impeller configurations. [START_REF] Petříček | Gas-liquid mass transfer rates in unbaffled tanks stirred by PBT: scale-up effects and pumping direction[END_REF] examined some of the relationships observed in the literature.

k L a = 0.023 P g V 0.44 U 0.47 g 1 + µ app µ c -1 . (3.31)
The proposed correlation obtained by Gabelle et al. successfully applied to reacting flows in stirred tanks [START_REF] Laakkonen | Modelling local gas-liquid mass transfer in agitated vessels[END_REF], is adopted in this work, due to its efficiency and adaptability to our case study, and also for its previous availability Fluent thermodynamics package.

Bubbles aspects and distribution in liquid sparged gas stirred reactors

Knowledge of bubble size in a gas-liquid reactor is critical because it allows for a better understanding of the gas dispersion process and reactor performance. Because of the complicated hydrodynamic environment, predicting the bubble size distribution in an aerated vessel is difficult. In reality, many amounts of turbulence can be encountered in a tank, particularly when the liquid has a complex rheology. The most widely used formula for estimating the mean bubble size inside a stirred tank was developed [START_REF] Ph Calderbank | Physical rate processes in industrial fermentations part i: The interfacial area in gas-liquid contacting with mechanical agitation[END_REF] d 32 ∝ σ 0.6 ε 2.5 ρ 0.6 α 0.5 g .

(3.32)

In fact, the Sauter mean diameter is directly related to the interfacial area a with the gas holdup:

a = 6 α g d 32 1 -α g . (3.33)
According to this formula, the Sauter mean diameter is correlated to the maximum diameter, which is determined by the ratio of turbulent shear stresses to interfacial forces: when the first surpass the second, the deformed bubble breakage. Equation (3.34) is thus based on turbulence theory and is true when no coalescence occurs; Calderbank's correlation is frequently seen in the following simplified version

d 32 ∝ P g V -0.4 . (3.34)
The exponent 0.4 follows from the assumptions that in turbulent flows bubble size is solely determined by the breakup probability [START_REF] Ss Alves | Bubble size in aerated stirred tanks[END_REF]. For coalescing systems, Alves and his team proposed the following empirical correlations

d 32 B = 0.0076 P g V -0.14 (3.35) d 32 T = 8.5 1 + 32.5 Q g D 2 P g V I -0.24 (3.36)
here the subscripts B and T indicate the mean diameter in the bulk and in the turbine region, respectively, whereas V I represents the volume swept by the impeller. Effect of gas loading was expressed in terms of gas flow rate Q g .

The interaction between the gas and liquid phases in the tank when gas is added can significantly affect the system's hydrodynamics. Three distinct regimes based on the bubble trajectories have been reported in the literature, and the resulting flow pattern depends on the operating conditions and the design of the chemical reactor. These regimes are: Flooding, loading, and fully recirculated. Figure 3.8 illustrates the change in flow regimes produced by modifying the impeller speed or gas flow rate. When the rotation speed is too slow for a given gas velocity, the liquid velocity is insufficient to push the gas, causing the impeller to flood. With such minimal gas phase dispersion within the STR, the gas creates ragged voids behind the impeller blades and effectively adopts helical trajectories around the shaft.

The top of the reactor started to fill with bubbles as the rotation speed rises, and the regime changes from flooding to loading. The gas phase is distributed evenly throughout the tank as a consequence of further increases in impeller speed. Gas bubbles can also be observed in the recirculation loops above and below the impeller disk when the impeller speed is very high, and extra gas recirculating vortexes can be detected close to the free surfaces.

On other hand, Nienow (1998) and al. found that gas tends to build up and generate perforated cavities in low-pressure areas behind the blades The existence of these cavities reduces the power generated, resulting in the impeller providing less energy to the liquid. The lower the power drawn, the deeper the cavities. Figure 3.9 depicts the size and shape of these cavities formed behind the blades of a Rushton turbine. The shape of the cavities may be predicted thanks to empirical correlations since cavity development is strongly related to the gas flow regime. Complex correlations link the gas flow number, Froude number, and stirred reactor shape [START_REF] Aw Nienow | Hydrodynamics of stirred bioreactors[END_REF].

Recent investigations on the breakup of single bubbles in stirred tanks have defined bubble breakup as a phenomenon caused primarily by turbulent fluctuations of momentum in the liquid phase; these are described as eddy-bubble interactions (li2; sha). Depending on the tur- [START_REF] Andersson | Modeling the breakup of fluid particles in turbulent flows[END_REF] 3.5. Turbulence-Chemistry Interactions (TCI) for multiphase react-

ing flows

The turbulent mixing process is extremely complex, owing to the fact that it includes a large variety of time and length scales. The mixing process begins at the macroscale level and descends down to the micro scale level. The macroscale level is essential because it serves as the origin of the mixing process. Chemical reactions, on the other hand, may occur only when the reacting species come into contact at the molecular, or microscale, level. The main challenges in modeling turbulent reacting flows are handling the mean rate of reaction and adequate representation of the chemistry in the model [START_REF] Fournier | A new parallel competing reaction system for assessing micromixing efficiency-experimental approach[END_REF].

Because the computational grid used for CFD calculations does not, and cannot, extend down to such a fine scale for applications of chemical engineering interest, it is impossible to include a detailed description of all relevant turbulence scales into a CFD code, particularly those that occur at the sub-grid or molecular level. A computational grid used to simulate a real physical system with a grid fine enough to capture molecule size phenomena would need computer resources well beyond the capabilities of today's computers.

In CFD modeling discipline, modeling the chemical turbulent kinetic is mature enough for combustion studies, with a very limited use in industrial manufacturing processes, since its classified as the advanced CFD modeling option (Manual, 2009). the diagram in figure 3.11 illustrates a classification of the developed Turbulence Chemistry Interactions models.

Following this diagram, and as previously mentioned the DAP preneutralization chemical reaction is a fast very reaction, the Eddy dissipation model based on the species transport in a multiphase flows in adopted. reactants will react as soon as they come into contact on a molecular level. Assuming that the EDM has a fast chemistry limit, the rate at which reactions occur is then controlled by the rate at which turbulent eddies carrying the reactants interactions. In these other expression, a turbulent mixing time primarily controls the mean response rate. This mixing time is calculated on a dimensional basis from the integral length scales using the turbulence model parameters that represent the energy cascade mechanism in turbulent flows.

This model is described based on the conservation equations, and for each phase k, it predicts the local mass fraction of each species, Y i k , through the solution of a convection-diffusion equation for the i th species. The equation is written then in the following form

∂ ∂t ρ q α q Y q i + ∇ • ρ q α q ⃗ v q Y q i = -∇ • α q ⃗ J q i + α q R q i + α q S q i + n p=1 ṁp i q j -ṁq j p i + R (3.37)
where R i q is the net rate of production of homogeneous species i by chemical reaction for phase q, ṁq j p i is the mass transfer source between species i and j from phase q to p, and R is the heterogeneous reaction rate. In addition, α q is the volume fraction for phase q and S i q is the rate of creation by addition from the dispersed phase plus any user-defined sources. The net rate of production of species i due to reaction r, R i,r , is given by the smaller (i.e., limiting value) of the two expressions below

R i,r = ν ′ i,r M w,i Aρ ϵ k min R Y R ν ′ R,r M w,R (3.38) R i,r = ν ′ i,r M w,i ABρ ϵ k P Y P N j ν ′′ j,r M w,j (3.39)
where Y P is the mass fraction of any product species, P , Y R is the mass fraction of a particular reactant, R. A is an empirical constant equal to 4.0, and B is an empirical constant equal to 0.5

Arrhenius rate, and reactions are mixing-limited (Han et al., 2020).

Hydrodynamics of stirred sparged gas liquid reactors

The overall characterization of a mixing systems is often performed by analyzing this main parameters and flow characteristics: the Reynolds number Re, power consumption P , power number N p , pumping capacity Q p , the shear rate γ and vorticity, mixing time Θ t , and vorticity.

In this section, these features will be discussed one by one.

Reynolds number Re & Froud number F r

The Reynolds number represents the ratio between inertial forces and forces viscous in a system (Equ. 3.40), it characterizes whether a flow is turbulent or laminar:

Re = ρ.N .D 2 µ . (3.40)
Thus, it is possible to make the distinction between the three flow regimes present in a mixing system [START_REF] Ottino | Mixing, chaotic advection, and turbulence[END_REF]):

• Laminar flow: Re<10, in this regime the viscous forces are predominant over the inertial forces. Based on Ottino's model, this regime is described in three stages: stretching, folding and breaking a phase;

• Transient regime: 10 <Re <1000; this scheme is, as its name suggests, the transition between the laminar regime and the turbulent regime.

• Turbulent regime: Re> 1000, in this regime the inertial forces predominate over the viscous forces. Moreover, this regime induces a random and asymmetric flow in the mixing system.

It is worth noting that transition between regimes in the gas liquid mixing system is also possible based on Froude number Fr, in the following formula:

Fr = N 2 .D g . (3.41)
The dimensionless number of Froude characterizes the importance of gravity in a mixture by representing the ratio between the forces related to the impeller speed and the force of gravity, with a specific gas flow rate. As the majority of hydrodynamics flow parameters, it depends directly to the impeller type, as illustrated in Figure 3.12 that represents the flow map for regimes classification within a mixing system that contains a Rushton turbine impeller [START_REF] Mmcg Warmoeskerken | Flooding of disc turbines in gas-liquid dispersions: A new description of the phenomenon[END_REF].

Figure 3.12: Flow map presentation for a Rushton turbine

Power consumption P

The power consumption by the mixing systems, otherwise called the "dissipated power" is expressed by the dissipated energy. It is a crucial characteristic that allows to compare the energy consumption.

P = 2πN τ, (3.42)
where, N is the impeller rotational speed, τ is the torque exerted on all blades was computed from the total momentum vector, which is computed by summing the cross products of the pressure and viscous forces vectors for each facet on the impeller with the moment vector.

For stirred sparged reactors, there are closely relationship between the gas-impeller design and the power demand. Figure 3.13 shows the dependence of the relative power demand (RPD) to the gas flow number for a STR equipped with a PBT impeller; as more gas is delivered to the impeller (or the impeller's speed decreases), there is a greater tendency for gas to accumulate in the low-pressure regions behind the blades, creating ventilated "cavities" (Warmoeskerken and Smith, 1986). When they are large, they can lead to a significant decrease in the impeller's power number (due to their blocking of the liquid discharge from the impeller) Ranganathan and Savithri (2010).

Power number N p

One of the most important dimensionless numbers for mechanical design, hydrodynamics evaluation of STRs and comparing the performance of different impellers is the "power number" (N p ), given by :

N p = P ρN 3 D 5 , (3.43)
where, P is the impeller power consumption (W), N is the impeller speed, ρ the density of the fluid, and D the impeller diameter [START_REF] Akiti | Turbulent mixing and chemical reaction in baffled stirred tank reactors: a comparison between experiments and a novel micromixing-based computational fluid dynamics model[END_REF]. Therefore, it is possible to establish a relationship between the power number and the Reynolds number in order to quantify the power consumption for a given mixing system and flow regime Fig. 3.14. Generally the N p is constant in turbulent regime. Pumping capacity is also an important factor when analyzing a mixing system. It allows to calculate the agitator's ability to produce a flow in the agitated reactor. It denotes the volumetric fluid flow through the volume swept by the system agitator. There are, however, many ways for evaluating the pumping created, depending on the shape of the agitators. The fundamental method involves constructing a box around the agitator and analyzing incoming and outgoing flows [START_REF] Guida | Pept measurements of solid-liquid flow field and spatial phase distribution in concentrated monodisperse stirred suspensions[END_REF]. Since the conservation of the mass must be respected inside this box (Incoming flow = outflow), it is then possible to evaluate the pumping rates generated by turbine type agitators. However, this technique does not work with complex geometry agitators.

To answer this problem, Devals et al. [START_REF] Devals | Cfd analysis of several design parameters affecting the performance of the maxblend impeller[END_REF] presented in the following formulation

Q + z (z) = 1 H H 0 ( A v + z dA)dz. ( 3
N q = Q N d 3 , (3.45)
In the literature, results for N q are less common than for the power number. (Dickey and JC, 1976) suggested a theoretical work in accordance with Nagata's (1975) experimental results for baffled and unbaffled vessels agitated by an 8-blade paddle impeller (Fig. 3.15). N q =0.75 was determined using LDV measurements in baffled tanks agitated by a Rushton turbine in totally turbulent flow. [START_REF] Yapıcı | Experimental determination of the optimum performance of ejector refrigeration system depending on ejector area ratio[END_REF] quantitatively replicated this results. Theoretical analysis proves that Nq must approach toward 0 for any shape since secondary velocities in creeping flow are insignificant. The current CFD findings and the numerical simulations of [START_REF] Takahashi | Chaotic mixing created by object inserted in a vessel agitated by an impeller[END_REF] proposal for a twin shaft mixing vessel both supported this behavior.

Figure 3.15: Pumping number for a baffled and unbaffled chemical reactor agitated by a RT [START_REF] Kresta | Advances in industrial mixing: a companion to the handbook of industrial mixing[END_REF].

Shear rate

Another important parameter in the characterization of a mixing system is the shear rate generated by the agitator. The shear rate is defined as a velocity gradient (Eq. 3.46)

γ = ∂v ∂y (3.46)
Indeed, the shear generated by the agitator is responsible for several phenomena often present in the mixture such as gas bubble dispersion, stretching and breaking phases and many other phenomena. It is therefore important to consider this phenomenon in the evaluation of a mixing system. However, the shear is not constant in a mixing tank: it is maximum at the end of the blade of the agitator and minimal at the wall of the tank. The non-dimensional shear count N s , defined by Equation (3.47)

N s = || γ|| N (3.47)
3.6.6. Mixing time Θ t

The dimensionless mixing time parameter is another global characteristic of a stirred reactors that has been traditionally measured by experiment, but which is now potentially accessible by CFD simulation. In experimental studies, the mixing is quantified by the time required for the concentration fluctuations of tracer material to have decayed to within 5% of the final wellmixed concentration (95% mixing time, θ 95 ) [START_REF] Devals | Cfd analysis of several design parameters affecting the performance of the maxblend impeller[END_REF]. Its expressed by

N θ 95 = f (Re) (3.48)
At high Reynolds numbers the dimensionless mixing time is a constant. There have been many attempts to calculate this mixing time, by solving the species transport equation, to predict the transient distribution of tracer that has been added to a developed flow inside a stirred tank.

Vorticity

The free surface in fluid mechanics is the boundary between a fluid and a gas, which is frequently present in liquid form. Through mixing tanks, a variety of subjects, such hydrology, the research of tsunamis, and the study of dams, approach the study of the free surface. When a fluid stops being at rest, the free surface of the system deforms. This deformation might be stationary or non-stationary, depending on the disturbance the fluid causes. The deformation of the originally flat free surface at rest is a special case of which the vortex is one. The free vortex is present when the angular momentum of the velocity is uniform throughout the fluid. The tangential velocity varies inversely with the distance from the center of rotation [START_REF] Aref | Stirring by chaotic advection[END_REF].

That is, the fluid away from the center of the vortex makes a revolution more slowly than the fluid near the center. Thus, the velocity of the fluid increases as it approaches the center of the vortex. The forced vortex corresponds to a solid rotation of the fluid volume. Not having shear, the rotational velocities acting on the fluid is identical everywhere. Thus, the tangential velocity takes the form given by the following equation and the free surface of the fluid adopts a parabolic form [START_REF] Aref | Stirring by chaotic advection[END_REF].

v θ = Γ 2.π.r (
v θ = ω.r (3.52)
In the case of agitated systems, developed vortices adopt a composite shape: a forced vortex form near the center and a free vortex near the tank walls. For our research work we are looking to study the vortex formation in the presence and absence of baffles inside the preneutralizer defined as a stirred tank with a non standard geometry.

Conclusion

In the fertilizer phosphate industry, the preneutralizer stirred reactor is an essential operation unit. That chapter presents the different aspects concerning the Stirred tank reactors, such as the geometric properties of agitated systems and their impact on the mixing efficiency, and also the hydrodynamic parameters. The CFD approach is also highlighted in this chapter including the mathematical models, the multiphase modeling approaches, the conservation of mass, momentum, and energy, as well as the turbulence effect. We presented the finite volume approach to solve these equations numerically. The modeling of the impeller rotation, the turbulence chemistry interactions modeling, and the the main hydrodynamics evaluation parameters are investigated. In the next chapter, theoretical background of CFD modeling in rotary drums will be investigated.

Introduction

Granulation is a complex process of controlling the desired properties of the particles by managing the process design and formulation composition. Regarding the required particle characteristics, the particle size, porosity, and moisture content are some of the most often used required particle characteristics.Its defined as the conversion of small, movable aggregates (particles) into big, semi-permanent aggregates (granules) that still allow the differentiation of the original aggregates. The two main types of granulation techniques are dry granulation and wet granulation, respectively. In dry granulation techniques, fluids do not enhance the granulation process. Wet granulation, as opposed to dry granulation, depends on a slurry (liquid binder or melt) to keep the granulation system going. In a granulator, the slurry is sprayed or poured onto the particles or onto the moving surface while the particles (fine granules) are being churned.

During the wet granulation procedures, the material is typically crushed and turned [START_REF] Adetayo | The effect of process parameters on drum granulation of fertilizers with broad size distributors[END_REF][START_REF] Litster | The science and engineering of granulation processes[END_REF].

The main objective of the industrial processing of phosphate minerals is the production of granular fertilizers, that are very soluble and directly assimilable by plants ysing the wet granulation process [START_REF] Dittmar | Ullmann's Encyclopedia of Industrial Chemistry[END_REF]. The fertilizer granulation essentially aims at forming agglomerates of a well-defined size that is generally larger than the initial particle size. The process also targets specific physico-chemical properties for the phosphate fertilizer grains, in occurrence, particle size distribution, density, porosity, solidity, and residual water content [START_REF] Artyukhov | Hydrodynamics of gas flow in small-sized vortex granulators in the production of Nitrogen fertilizers[END_REF]. The wet granulation process is based on different technologies with different designs and mechanisms depending on the physico-chemical properties of the input product [START_REF] Saleh | Mise en oeuvre des poudres: Techniques de granulation humide et liants[END_REF]; [START_REF] Ntamo | Industry 4.0 in action: Digitalisation of a continuous process manufacturing for formulated products[END_REF]. fertilizers industry, they are still challenging equipment that require a complex adaptation of the operating conditions to the quality of the raw material feeding the process [START_REF] Yang | Numerical simulation of particle dynamics in different flow regimes in a rotating drum[END_REF].

Industrial granulation technologies

Depending on the operating parameters, the physical properties of materials, and the geometrical designs, rotating drums can show different flow regimes for the particles dynamics Govender (2016). The flow regime significantly affects the quality of material processing in terms of mixing, grinding, size enlargement, etc. Experimental studies conducted by [START_REF] Arntz | Granular mixing and segregation in a horizontal rotating drum: A simulation study on the impact of rotational speed and fill level[END_REF] showed that rotating drums can exhibit seven flow regimes : Sliding, surging, slumping, rolling, cascading, cataracting and centrifuging [START_REF] Pathmathas | Granular flow modelling of rotating drum flows using positron emission particle tracking[END_REF]. Each regime is characterized by its specific hydrodynamics Santos et al. (2016b).

Complex phenomena such avalanching, segregation, mixing, drying, aggregation, comminution (abrasion, attrition, and fracture), and convection occur during the movement of particles in spinning drums. In industries, mixing, rotary kilns, and comminution all use the flow of granular material in revolving drums. While mixers and kilns typically require an active flowing layer to mediate the mixing or heat transfer requirements, the latter application (tumbling mills) typically requires high frequency rotations of dense media in drums fitted with radial baffles to facilitate the strain rate and cataracting demands of the process.

Granulation mechanisms

Granulation is a complex process with several competing physical phenomena occurring in the granulation vessel, which lead ultimately to the formation of the granules [START_REF] Litster | The science and engineering of granulation processes[END_REF]). We will divide these phenomena into three groups of rate processes (figures 4.3-4.5):

• Wetting, nucleation and binder distribution The first stage of granulation is the addition and distribution of the binder to give nuclei granules. This occurs within the liquid spray zone in the granulator. The liquid binder is usually sprayed onto the moving powder bed. Ideally, each drop will imbibe into the powder bed, engulfing particles to form a single granule nucleus. If the drop does not easily wet the powder, or the rate of imbibition is slow, large wet agglomerates will format the powder surface.

Sometimes shear forces within the powder bed are large enough to break up these clumps of wet material to further distribute the liquid.

Nucleation gives a distribution of loosely packed granule nuclei. Under good wetting conditions, the distribution will be narrow and correlate closely with the liquid drop size distribution.

Poor wetting and nucleation causes a broad nuclei size distribution and often poor final granule properties. The nucleation process is strongly influenced by the characteristics of the liquid spray, the flux of powder through the spray zone and the wetting properties of the formulation.

Granule nuclei will consolidate through collisions with other granules, and granulator. The extent of consolidation depends on the intensity of agitation in the granulator and resistance of the granule to deformation. Granules made from fine powders or with viscous liquid binders resist deformation on impact and consolidate slowly. Granules made from coarse powders deform a lot on impact and quickly reach their minimum porosity. Granule consolidation controls the final granule porosity which influences many other granule properties.

When two granules collide they may stick together to form a single large granule. This is growth by coalescence. For successful coalescence (a) the energy of impact must be absorbed during collision so that the granules do not rebound; and (b) a strong bond must form at the contact between the colliding granules. The presence of liquid at the surface of the granule is important for growth by coalescence and coalescence rate is very sensitive to liquid content. Consolidation and coalescence are closely related and are considered together in this book. Granules formed through a combination of nucleation and growth processes are eventually dried. The dry granules may break to give several fragments or attrit to give a fine powder. These processes are controlled by the fracture mechanics of the granule. Granules resist attrition if their fracture toughness is high (related to formulation properties) and the number and size of flaws is small (related to granule porosity). These rate processes are at the heart of the granulation process and determine the final granule properties. It is important to know which rate processes dominate in any particular application. For example, if nucleation is the dominant rate process, adjusting the liquid spray nozzles will have a profound effect on granule size distribution. If coalescence rates are high, granule size will be very sensitive to liquid content and granule residence time.

We therefore spend a lot of time examining each of these processes in detail. We will:

• examine the underlying physics behind each rate process

• define the controlling formulation properties and process parameters for each rate process

• use regime maps to establish the operating regime for the granulator

• provide quantitative relationships to predict the effect of changing operating parameters and formulation properties

• provide troubleshooting and scale up guidelines

Wetting and nucleation

The process of granulation starts with wetting and nucleation (Figure 4.3). The liquid spray droplets and the granules interact in the spraying zone. The liquid is being uniformly spread across the granules. The first aggregates are created as a result of this step. Granule formation is significantly influenced by wetting and nucleation [START_REF] Simon M Iveson | Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review[END_REF][START_REF] It Cameron | Process systems modelling and applications in granulation: A review[END_REF]. However, this mechanism is hardly ever recognized and distinguished from the other two granulation mechanisms due to an inadequate understanding of the variables that govern the wetting and nucleation processes.

Growth and consolidation

Particle layering and particle aggregation are the two main processes that allow particles to grow (Figure 4.4). When a liquid phase is continuously coated onto a particle surface, continuous par- ticle growth also known as particle layering-occurs. The surface of the particle forms a "onion skin" once the liquid phase, which can take the form of a melt, solution, or slurry, solidifies.

The particle's volume and mass consequently grow. However, The system's particle count stays constant. This process results in the production of compact and hard granules. A discrete particle growth method is known as particle agglomeration. Such particle growth behavior is carried on by particle collisions and the adhesion of the successful collisions. Binary particle agglomeration happens when two particles successfully contact and agglomerate together to produce a larger, composite particle (Figure 4.4). Particle agglomeration can occur both early in the granulation process (simultaneously with particle wetting and nucleation) and subsequently in the process, after the liquid has already been injected to the granulator [START_REF] Simon M Iveson | Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review[END_REF][START_REF] Mort | Scale-up of agglomeration processes using transformations[END_REF]. The number of particles is decreased as a result of binary agglomeration since two particles die and one new particle is born.

Attrition and breakage

Wet particle breakage is a discrete process that affects the system's particle population (Figure 4.5). In high shear equipment, particularly in high intensity mixer granulation, breakage effects are crucial. On the other hand, the particle size is substantially less affected by the particle breaking phenomenon in rotary drum granulators than by the growth process. In light of this, limited investigations have been done on particle breakage in rotary drums. Breakdown of dry particles is called as particle attrition. The particle layering mechanism's impact is reversed by the attrition rate, which is a negative rate. The particle attrition rate is strongly influenced by the particle velocities inside the granulator and is proportional to the bed conditions. In the granulation operations, dust generation is a common result of particle attrition and is thus extremely undesired [START_REF] Askarishahi | Mechanistic modelling of fluid bed granulation, part i: Agglomeration in pilot scale process[END_REF]. However, granulators that employ high granule velocities need to consider this mechanism.

Numerical modeling approaches for multiphase flows with granular phase

Granulation processes need to be mathematically modelled in order to operate well. Both offline optimization and online control and monitoring employ mathematical models. Depending on the goals of the study, several modeling strategies for granulation processes could be chosen.

Nowadays, the increase of computing capabilities makes numerical simulation quite affordable to study the particle dynamics in rotating drums. Several studies have proven the CFD efficiency for granular flows assessment and analysis. For these simulations, the Euler-Euler The Euler-Euler approach considers the fluid and solid phases as interpenetrating continua, coupled through interaction terms. The solid phase properties are obtained using the kinetic theory of granular flow. The relative computational efficiency of this approach with respect to the other ones cited above, as well as the fast development in computational resources allow it to replace, at least partially, the experimental studies, while giving a way to analyze phenomena that are hardly studied experimentally. Yet, great attention should be given to the interaction models.

In Particular, the interphase momentum exchange models are a key component for an accurate

Population balance modeling (BPM)

Population balance modeling (BPM) is frequently used in the following applications: I process design to predict particle properties, such as predicting the produced particles size distributions (PSD); (ii) process control and optimization; and (iii) sensitivity analysis to investigate potential changes in product characteristics quality brought on by a change in operating conditions. The usefulness of PBM approach in granulation processes is nevertheless constrained by inadequate comprehension of the rate processes (quantitative representations of granulation mechanisms)

and difficult solution techniques resulting from the PBMs, particularly for agglomeration issues.

Th concept of population balance method (PBM) was defined by Ramkishna as an equation in the number density, which means the number balance of particles for a particular state [START_REF] Ramkrishna | Population balances: Theory and applications to particulate systems in engineering[END_REF]. In the PBM model equations, particle properties are denoted as an internal coordinate, while particle location (spatial position) is denoted as an external coordinate. Some of the particle properties (internal) of a special interest are: size, shape, density, moisture content, gas content, porosity, composition, age, etc. Mathematical representation of the number-based PBM with one internal (particle diameter), and one external (spatial position in the granulator) coordinate is given by

∂n(d, z, t) ∂t = - ∂ ∂d [Gn(d, z, t)] + B a (d, z, t) -D a (d, z, t) + B b (d, z, t) -D b (d, z, t) - ∂ ∂z [vn(d, z, t)], (4.2)
where n(d, z, t) is the number density function, d is the particle diameter, z is the spatial position, and t is the time. The first term on the rhs. represents the particle growth due to layering at a growth rate G; the second and the third terms on the rhs. Respectively, for particle birth (B a ) and death (D a ) due to agglomeration. The particle birth (B b ) and death (D b ) rates due to breakage are represented by the fourth and the fifth terms on the rhs, respectively. The final term on the rhs represents a continuous process and gives the flow of particles with the velocity v through the granulator. The layering term G, as previously mentioned, is expressed in the following.

The creation of a coating of the slurry on the particle provides as the stimulus for layering, a continuous process (differential growth). A particle size independent linear growth rate is frequently assumed in PBM models. Growth rate independent of particle size.

Representation of the agglomeration terms, B a and D a : particle agglomeration is a discrete event that results in reduced particle number in the system. Quantitative representation of the agglomeration process is very challenging. One of the most widely used agglomeration process models for granulator-scale modelling was introduced by Hulburt and Katz in [START_REF] Ramkrishna | Population balances: Theory and applications to particulate systems in engineering[END_REF]. Formulation of the agglomeration process assumes a binary particle agglomeration process, i.e., two particles colliding with each other. During binary particle agglomeration, two particles die, and one particle is born [START_REF] Ramkrishna | Population balances: Theory and applications to particulate systems in engineering[END_REF]. The birth of particles due to binary agglomeration is represented by

B a = 1 2 v 0 β(t, v -w, w)n(t, v -w)n(t, w)dw. (4.3)
The above equation shows the birth of particles of volume v as a result of agglomeration of particles of volumes vw and w. Merging of particles of volume v with any other particle, i.e., the death of particles is represented by the equation

D a = -n(t, v) ∞ 0 β(t, v, w)n(t, w)dw. (4.4)
In the previous equations, β is the agglomeration rate, also called the agglomeration kernel. The agglomeration kernel defines the collision frequency of the two colliding Representation of the breakage terms, B b and D b , particle breakage in rotary drum granulators is not significant [START_REF] Ramkrishna | Population balances: Theory and applications to particulate systems in engineering[END_REF], therefore, only a few studies regarding wet breakage in rotary drums granulators have been investigated.

Discrete Element Modeling (DEM)

DEM is a Lagrangian method for modeling individual trajectories of every particle in a granular system. Particle trajectories are calculated by solving Newton's equations of motion, in which contact forces and cohesive forces are considered. The movement of an arbitrary particle i, which is in contact with a set of j = 1, . . . , k adjacent particles, is described by

m i ẍi = k j=1 F c ij + F l ij + m i g, (4.5
)

I i θi = k j=1 h ij × F c ij + F l ij + q ij , (4.6)
with mass m i , position x i , contact force F c ij , cohesive force F l ij , gravity g, mass moment of inertia I i , angular position θ i , distance-to-contact point vector h ij and additional torque q ij (due to rolling). Note that variables printed in bold represent vectors. A linear spring-dashpot model was used to compute the total frictional force F c ij between two colliding particles i and j. F c ij is nonzero only if there is a positive absolute normal overlap δ n between the particles, which is defined as

δ n = r i + r j -x i -x j • n ij , ( 4.7)
where r is the particle radius. When δ n > 0, F c ij is computed as

F c ij = k n δ n + γ n δn + k t δ t + γ t δt , (4.8)
with elastic constant k, overlap δ and visco-elastic damping constant γ. The subscripts n and t denote normal and tangential directions respectively. Furthermore, δ t is truncated to fulfill

k t δ t + γ t δt ≤ µ F k n δ n + γ n δn , ( 4.9)
with coefficient of friction µ F , to ensure the satisfaction of Coulomb's law of attrition. The Hertz-Mindlin model, based on Hertz theory in normal direction and the Mindlin no-slip improved model in tangential direction,is genreraly employed to calculate the elastic and damping constants. k n is calculated by

k n = 4 3        1 -ν 2 i E i + 1 -ν 2 j E j        -1 R * δ n , (4.10)
with modulus of elasticity E, Poisson's ratio ν and effective radius R * .γ n is computed as

γ n = - 2 √ 15 3 ln(e) ln 2 (e) + π 2        1 ρ 1 πr 3 i + 1 ρ 2 πr 3 j        -1 k n (4.11)
with coefficient of restitution (CoR) e and particle density ρ. As for the tangential components, k t is defined as

k t = 8 2 -ν i G i + 2 -ν j G j -1 R * δ n (4.12)
DEM approach allows to consider the liquid phase with a minimal fraction, from Easo et al.

proposition [START_REF] Kloss | Models, algorithms and validation for opensource dem and cfd-dem[END_REF], based on the definition of cohesive forces between two particles i and j. It is assumed that a surface liquid film exists on each particle, allowing for a liquid bridge to form upon contact with another particle. For the purpose of the cohesion model, particles are considered to be in contact only if δ n satisfies:

-0.1 r i + r j ≤ δ n ≤ -0.01 r i + r j .

(4.13)

The bridge ruptures when x ix j exceeds the rupture distance d 0 , which is given by

d 0 =         1 2 + α i + α j 4         V 1/3 bond , (4.14)
with contact angle α and liquid bond volume V bond . It is assumed that V bond is 5% of the combined liquid film volume, and that distributes evenly over the two particles after bridge rupture. The total cohesive force F l ij between the particles is defined as with liquid viscosity µ. In the case of particle-wall contact, F cap and F vis are equal to zero.

F l ij = F cap + F vis , ( 4

CFD modeling of granulation (CFD + KTGF)

Both Eulerian and Lagrangian methods may be used to describe multiphase granular flows in rotating drums computationally. The Discrete Element Method (DEM), which belongs to the Lagrangian family, works by simulating collisions based on interaction rules and then predicting the trajectories of each individual particle. Many publications have utilized this method to characterize the angle of repose for various regimes and to describe mixing processes in order to describe the bed behavior in rotating drums.

The relevant works are based on the Eulerian-Eulerian formulation of the balance equations coupled with the kinetic theory for dense granular flows (KTGF) inside rotating drums Santos et al. (2016b); [START_REF] Santomaso | Investigation of the granular behaviour in a rotating drum operated over a wide range of rotational speed[END_REF]; [START_REF] Khalilitehrani | Modeling dilute and dense granular flows in a high shear granulator[END_REF]. To predict such complex flows, a lot of effort has been put into improving numerical methods such Computational Fluid Dynamics (CFD) tools. However, it has been shown that multiphase flow modeling systems with one or more particle phases are the most complicated and difficult. It is crucial to select a numerical technique that can account for particle-wall and particle-particle interactions in three dimensions and across any particle size distribution in order to properly predict the behavior of solids. The prediction of gas-solids flows is achievable using a variety of CFD models.

The Euler-Euler CFD modeling approach have been presented with details in chapter 3. This section described the adapted Euler-Euler model to a multiphase flow with a granular phase.

Neglecting the effect of compressibility, the continuity equations for gas-solid mass conservation are represented by the following equations

∇ • α g v g ϱ g = 0, ∇ • (α s v s ϱ s ) = 0, (4.18)
where, the subscripts g and s stand for the gas phase (air) and the solid phase (particles), respectively. α, ϱ, and v stand for the volume fraction, density and velocity, respectively, with α g + α s = 1. The conservation of the gas and solid momentum is given by the following equations

∂ ∂t α g ϱ g v g + ∇ • α g ϱ g v g v g = -α g ∇p+ ∇ • τ g + α g ϱ g g + α g ϱ g F ext g + F lift g + F vm g + K gs v s -v g , (4.19) ∂ ∂t (α s ϱ s v s ) + ∇ • (α s ϱ s v s v s ) = -α s ∇p + ∇p s + ∇ • τ s + α s ϱ s g + α s ϱ s F ext s + F lift s + F vm s + K gs v g -v s , (4.20)
where F ext q (q stand either for g or s), F lift q , and F vm q are the external body force, the lift force and the virtual mass force, respectively. τ q are the phase stress tensors given by :

τ g = α g µ g ∇v g + ∇v ⊺ g τ s = α s µ s ∇v s + ∇v ⊺ s + α s λ s -2 3 µ s ∇ • v s I (4.21)
K sg is the momentum drag exchange coefficient between the phases. In the current investigation, we assume that the system is isothermal, and only drag and gravity forces are considered.

Indeed, many research studies Santos et al. (2016a); [START_REF] Delele | Studying the solids and fluid flow behavior in rotary drums based on a multiphase CFD model[END_REF] concluded that the effects of the lift and virtual Mass forces are negligible compared to the drag force impact. We confirmed this with three simulation cases. The momentum exchange coefficient K sg between solid phase s and fluid phase g can be written in the following general form

K sg = α s ϱ s f d ϖ s ,
where f d is the drag function and ϖ s is the particulate relaxation time defined as

ϖ s = ϱ s d 2 s 18µ g ,
where d s is the diameter of particles in the solid phase and µ g is the air viscosity. The definition of the drag function f d includes a drag coefficient C d that is based on the relative Reynolds number Re r defined as The first one was developed by [START_REF] Wen | Chemical engineering progress symposium series[END_REF] and is given by

Re r = ϱ f v g -v b d s µ g . ( 4
K sg = 3C d α s α g ϱ g v s -v g 4d s α -2.65 g (4.23)
where

C d = 24 α g Re s 1 + 0.15 α g Re s 0.687 . ( 4

.24)

The second one was introduced by Gidaspow (1994) to correct [START_REF] Wen | Chemical engineering progress symposium series[END_REF] 's model following the expression

K sg =          3C d α s α g ϱ g |vs-vg| 4d s α -2.65 g for α g > 0.8, 150α s (1-α g )µ g α g d 2 s + 1.75ϱ g α s d s v s -v g for α g ⩽ 0.8. (4.25)
The third considered model is due to [START_REF] Syamlal | The derivation of a drag coefficient formula from velocity-voidage correlations[END_REF] and it can be expressed as: The Kinetic Theory for granular flows is a theory that extends the ideas and tools applied successfully for the kinetic theory of gases to flows characterized by dispersed particles. The dispersed particles are treated statistically, characterized by a frequency distribution of velocities and collisions, and the fluctuating velocities are then related to the shear gradient. The final outcome of the application of this model to the granular flow is a transport equation that describes the behaviour of a "granular temperature" Θ associated to the random movement of particles whose diffusion coefficients are given by (several different) empirical correlations. To model the particle-particle interactions, the KTGF is used to define the stress term τ [START_REF] Khalilitehrani | Modeling dilute and dense granular flows in a high shear granulator[END_REF]. This theory is based on an analogy with the classical kinetic gas theory. For instance, a term θ s , representing the granular temperature, also called "pseudo-thermal", accounts for the kinetic energy associated with the particle velocity fluctuations [START_REF] Chao | Investigation of the particle-particle drag in a dense binary fluidized bed[END_REF]. The transport equation for the granular temperature is expressed as following

K sg = 3C d α s α g ϱ g 4V 2 r,s d s Re s V r,s v s -v g , ( 4
3 2 ∂ ∂t (ϱ s α s θ s ) + ∇ • (ϱ s α s v s θ s ) = (-p s I + τ s ) : ∇v s + ∇ • k θ s ∇θ s -γ θ s -3K sg Φ gs , , (4.30)
where (-p s I + τ s ) : ∇v s represents the resultant energy from the solid stress tensor, ∇ • (k θs ∇θ s ) is the diffusion of energy, Φ gs is the energy exchange between the fluid and solid phases, k θ s is the diffusion coefficient, and γ θs is the collisional dissipation energy given by Lun (1991):

γ θ s = 12 1 -e 2 ss α 2 s ϱ s g 0 θ 3/2 s d s √ π , ( 4.31)
and

θ s = 1 3 u 2 s,i , (4.32)
where e ss is the restitution coefficient for particle collisions, taken as 0.95, and the radial distribution function g 0 represents the probability of interparticle contact. Neglecting the convective and diffusive terms in Eq.(8.13) allows to reduce the transport equation for the granular temperature to the algebraic form

(-p s I + τ s ) : ∇v s -γ θ s -3K sg Φ gs = 0. (4.33)
To estimate the viscosity of the solid phase µ s , three terms are used:

µ s = µ fr s + µ col s + µ kin s , (4.34)
where µ fr s is the frictional viscosity, µ col s is the collisional viscosity, and µ kin s is the kinetic viscosity, which are modeled using detailed correlations that can be found in [START_REF] Schaeffer | Instability in the evolution equations describing incompressible granular flow[END_REF].

Hydrodynamics of granular flows in RDG: Regimes classification

Several experimental works and mathematical models were developed to predict the transitions between the different forms of transverse motion of rotating drums. Mellmann [START_REF] Mellmann | The transverse motion of solids in rotating cylinders-forms of motion and transition behavior[END_REF] identified this regimes classification of flow configurations for granular material in rotating drums depending on the filling level f, wall friction, and Froude number (Fr). According to Mellmann's classification there are three main regimes (Figure 4.6):

• Slipping: No granular flow occurs when rotational velocity, wall friction, and filling degree are extremely low. In the first, the material continually slips away from the wall, yet an inertial observer would not see any movement in the granular bed. Because wall friction is sufficient to "grip" the material in the latter case and to produce a periodic alternation between adherence and sliding from the wall, there is a bulk movement.

• Tumbling/Cascading; a "liquid-like" flow begins at intermediate levels of Fr. The slumping sub-motion, which is characterized by periodic cascades in the top of the bed that reduce the free-surface angle from its maximum value to its minimum value, can be seen as the rotational speed of the bed keeps increasing. Thereafter, the rigid rotation of the bed restores the higher angle, and the cycle repeats. Due to this, the almost flat free surface periodically oscillates between two limiting angles in the slumping regime. [START_REF] Mellmann | The transverse motion of solids in rotating cylinders-forms of motion and transition behavior[END_REF] • The rolling and cascading sub-regimes, both of which have an ongoing cascade, increase velocity much more. The shape of the free-surface, which is flat or S-shaped dependent on the motion, is the only distinction between rolling and cascading movements. This "liquid" regime is the most frequently seen in applications, as can be seen from the last row of figure 4.6, as it improves mixing.

• Cataracting: When the gravitational acceleration equals inertia, the particles are flung in a mass into th()e gas drum's empty space, initiating the cataracting sub-regime. The centrifuging process rotates the inside walls in a homogeneous layer of grains at even greater Fr.

Conclusion

Rotary drum granulator is a crucial component in phosphate fertilizer industry, since its the principal unit operation that govern the granules physico-chemical properties. That chapter presents a survey on the existing granulation technologies, especially for the wet granulation.

The granulation mechanisms from the wetting to the breakage are presented. Also, the used numerical models employed to the study of granulation processes are investigated, with a highlight on the PBM, DEM and CFD modeling approaches. A focus on the CFD KTGF coupling is presented, followed by the hydrodynamics characterization of the flow configurations. In the next parts, the results obtained from that research work will be investigated.

Part II

Preneutralization of Phosphoric Acid by Ammonia: Study, modeling and simulation using CFD Approach

Introduction

Various research studies [START_REF] Khamar | The use of the CFD for the hydrodynamic flow diagnostic and study in a Phosphoric Acid Reactor[END_REF] have been conducted in order to provide interesting information allowing to understand the complexity of the flow behavior in stirred tanks in phosphate fertilizer industry . This work uses CFD simulations of a preneutralizer chemical reactor characterized by a non-standard design in terms of the different characteristic dimension ratios of the vessel, and the shape. This stirred reactor is used in fertilizers manufacturing process as a preneutralization operation unit allowing to produce a chemical slurry with specific physicochemical characteristics and under special operating conditions [START_REF] Driss | 5-0028-0049A study of mixing structure in stirred tanks equipped with multiple four-blade Rushton impellers[END_REF]. The main goal of this chapter is to study the multiphase flow and mixing within the preneutralizer; and evaluate the effect of baffling on its hydrodynamic which should improve agitation performance and chemical species mixing.

Conceptual model of preneutralizer

CAD model and computational approach

Using a computational fluid dynamics (CFD) package, we had defined the computational domain equipped with four baffles spaced uniformly around the preneutralizer periphery as shown in Figure 5.1. The mixing system is a pitched blade turbine (PBT) with inclined blade impellers.

The CAD file is developed for the total reactor as presented in Figure 5.1, and in order to study the hydrodynamic flow behavior, we had activated only the part that will be full of the mixture (5.1.C). On the other hand, two computer aided designed (CAD) models are developed: Baffled and Unbaffled geometries. • A rotor defined by a rotating zone which surrounds the agitator;

• A stator as the stable and steady zone of the rest of the preneutralizer.

Hence, as explained in chapter 3, the Multiple Referential Frame (MRF) approach is based on the flow domain splitting in two parts. The first one is around the agitator that is simulated in a rotating frame with the stirrer and the equations are solved taking into consideration a Coriolis term. For the rest of the domain, the equations are solved in the fixed frame [START_REF] David | CFD modelling of stirred tanks: Numerical considerations[END_REF][START_REF] Hadane | Hydrodynamic study of a phosphate flotation cell by CFD approach[END_REF]. The MRF approach does not take into account the relative movement of a moving zone with respect to adjacent zones (which can be moving or stationary), the mesh remains motionless for calculation (Atiemo-Obeng et al., 2004). This is analogous to blocking the movement of the moving part in a specific position and the field of instantaneous flow around the rotor at this position is observed. This approach is an approximate one but it has the advantage of involving stationary calculations (the mesh is immobile) and therefore faster.

Domain meshing

A mesh was generated to discretize the domain into small control volumes, where the governing equations (continuity and momentum) were be approximated by numerical calculations. The mesh wash generated using ANSYS Mesher as presented in Figure 5.2.

Mathematical model and calculation setup

In the current study, we have considered that:

• The mixture behaves like an incompressible Newtonian viscous fluid,

• The system is isothermal (Temperature is constant throughout the sheath fluid of the preneutralizer reactor).

Operational inputs of the preneutralizer

Some parameters are crucial for calculating and simulating mixing in the preneutralizer, which we summarize in the following table 5.1.

Turbulence model

The aim of turbulence simulation is to predict the physical behavior of turbulent flow generated in a system using numerical methods [START_REF] Bakker | A computational model for the gas-liquid flow in stirred reactors[END_REF]. A modeling method used should ensure accuracy, simplicity, and computational efficiency . In order to study the flow behaviour and the mixing quality in the preneutralizer, the chosen model is based on the K-ε RNG, from Reynolds Averaged Navier Stokes Equations (RANS) two-equations turbulence models family, as described with details in chapter 3. This model is indeed adapted to industrial scale problems, and is a method where the smallest eddies are first resolved in the inertial range and then represented in terms of the next smallest eddies [START_REF] Louis | Flow and turbulence in a stirred tank[END_REF].

More details on setting the solver for simulating the preneutralizer are given in the following modeling capabilities, both with and without reactions. For this study we will use the transport species approach, to mention that we have seven equations to solve 1. Mass conservation equation:

∂ρ ∂t + ∇.(ρ. - → v ) = S m , (5.1)
the equation (5.1) is the general form of the mass conservation one and is valid for incompressible as well as compressible flows. The source S m is the mass added to the continuous phase.

2. Momentum Conservation Equations (Three dimensional system)

∂ρ. - → v ∂t + ∇.(ρ - → v - → v ) = -∇p + ∇( τ) + ρ - → g + - → F , (5.2)
With p is the static pressure, τ is the stress tensor, and ρ⃗ g and ⃗ F are the gravitational body force and external body forces.

For the turbulence model, as we had mention in section 3.3.4, the useful one is RNG k-ϵ, so we have two equations to solve :

3. Turbulent kinetic energy equation:

∂(ρk) ∂t + ∂(ρ.k.u i ) ∂x i = ∂ x j [(µ + µ t σ k ) ∂k ∂x j ] + P k -ρ.ϵ. (5.3) 4.
Turbulent dissipation equation:

∂(ρϵ) ∂t + ∂(ρ.ϵ.u i ) ∂x i = ∂ x j [(µ + µ t σ ϵ ) ∂ϵ ∂x j ] + C 1ϵ ϵ k P k -C * 2ϵ ρ ϵ 2 k , (5.4)
Where

C * 2ϵ = C 2ϵ + C ν η 3 (1- η η 0
) 1+βη 3 and η = Sk/ϵ , with S is the modulus of the mean rate of strain tensor.

Transport equations:

A prediction of the local mass fraction of each species (H 3 PO 4 and NH 3 : Y i , through the solution of a convection-diffusion equation for the i th species. This conservation equation takes the following general form:

∂ ∂t (ρY i ) + ∇.(ρ - → v Y i ) = -∇. - → J i + R i + S i (5.5)
Where R i is the net rate of production by chemical reaction and S i is the rate of creation by addition from the dispersed phase if it exist. For this stage of the study, we don't take into account the reaction between the defined species.So we have :

∂ ∂t (ρY H 3 P O 4 ) + ∇.(ρ - → v Y H 3 P O 4 ) = -∇. - → J H 3 P O 4 (5.6) ∂ ∂t (ρY N H 3 ) + ∇.(ρ - → v Y N H 3 ) = -∇. - → J N H 3 (5.7)
J i is the diffusion flux of species i, which arises due to concentration gradients, For this case we simulate a turbulent flow , so J i is defined by :

- → J i = -(ρ.D i,m + µ t Sc t )∇Y i (5.8)
With D i,m is the diffusion coefficient for species i in the mixture,µ t the turbulent viscosity, and Sc t is turbulent Schmidt number.

6. The seventh equation the Energy equation, for this study we suppose that we have an isothermal system , so this one doesn't make part of the system to solve.

Boundary conditions

To define a particular problem with a unique solution equation and transfer system, it is imperative to provide information on the flow at different boundaries of the domain. In this study, we have used the commercial solver ANSYS to solve the numerical model.

Thus, we have considered that:

• The liquid level surface (The top of the geometry) is modeled as a symmetry surface, which mean that the free surface is assumed to be flat, this is generally the good assumption for stirred vessels with at least three baffles;

• The agitator is divided in two parts, the first one is in the fixed reference frame and the second part is in the rotating reference frame;

• All the other boundary conditions are defined as standard walls.

Results analysis

The simulation is carried out for two cases, the first one is a baffled preneutralizer with the characteristics mentioned before, and the second case is for an unbaffled one, the reason behind is to evaluate the effect of baffles on the vortex formation within the preneutralizer. However, the vortex appearance in a phenomenon that impacts directly the hydrodynamics within stirred tanks, and for the reacting flows, the vortex inhibit the good mixing and restricts the best distribution of the chemical species along the STR.

Validation of flow regime

As explained in the section 3.3.5, we can make the difference between three regime flows: laminar, transient, and turbulent, according to the value of Reynolds number. For this study, R e = 38144695.05 which involved that the flow regime induced by the stirrer is turbulent. The choice of the model k-ε (RNG) is quite right.

Type of flow induced by agitation

The analysis of the shape of the currents lines developed within the preneutralizer fluid sheath in the two studied cases (Figure 5.3) shows that the fluid is explored from the blades of the agitator turbine towards the walls of the vessel and is divided into two parts, thus giving two loops of circulation which develop, one below. Consequently, it can be directly concluded that this stirrer ensures a radial flow of the mixture within the preneutralizer. 5. Baffles effects on the hydrodynamics within the preneutralizer using CFD approach specific cases depending to the geometry design modifications and the impeller rotational speed.

In the present case, the operating rotational speed is equal to 68 rpm, moderate to the domain dimensions. However, this rotational speed is chosen to ensure the required conditions for the chemical species interphase exchange. More than that, the non standard geometry design, with a tightening of the bottom diameter of the reactor, are behind promoting the radial flow near to the mixing zone.

Balance of the torques exerted on the agitator

For a speed of rotation of 68 rpm, the pairs resulting from the CFD are given in the following tables:

Power calculated using CFD

The general formula for the power is given by :

P = 2πN C (5.9)
With N; the number of rotations per minute (Hz), and C; the torque (N.m). For a speed of 68 rpm, the powers calculated using the CFD are given in the following table 5.4

Number of power N p

This parameter is dimensionless, and it's based on the hydraulic power calculated before, using the following formula (Equation 5.10)

N p = P ρN 3 D 5
(5.10)

The results obtained using the CFD post processing tools gives a value of N pu nbf = 0.2948 for the unbaffled preneutralizer and N pbf = 0,5559 for the baffled one. These values validate the previously published works, and proves the strong link between baffles presence and the hydrodynamical power within the flow. Four equal distance baffles allows to double the power and the power number, which enhances considerably the species dispersion in the system.

Pumping rate or Pumping capacity Q p

The pumping rate is the volume flow rate of fluid that passes through the agitator turbine.

Indeed, the stirrer rotating in the tank can be considered as a pump. This flow is also called the agitator pumping capacity noted Q p . From the results of the simulation Q punbf = 79.72 m 3 /s , and for the baffled preneutralizer the flow is Q punbf = 85.97 m 3 /s.

Number of pumping N q

It is the dimensionless formula of the pumping rate of an agitator. it depends essentially on the type of impellers and the hydrodynamic flow regime. It is given by the following expression:

N q = Q p N D 3 (5.11)
With Q p is the pumping flow (m 3 /s) and N is the rotation speed of the stirrer (rad/s), D is the turbine diameter (m). The results are N qunb = 0.51 and for the baffled one N qbf = 0.54 . We conclude the negligible effect of baffles on the pumping capacity and pumping flow rate, which explain its rate use to the hydrodynamical evaluation. Otherwise, this parameter is powerful for this analysis categories when the impeller type is modified. For our case, as we used the original PBT impeller, the variation for the two simulated cases is less significant.

Vortex detection

The results presented in Figure 5.4 shows the vortex form in the preneutralizer using the Qcriterion parameter [START_REF] Guyot | A Vortex modeling with 3D CFD[END_REF], that define the vortex as a connected fluid region with a positive second invariant. This criterion also adds a secondary condition on the pressure, requiring it to be lower than ambient pressure in the vortex. The Q represents the local balance between shear strain rate and vorticity magnitude, defining vortices as areas where the vorticity magnitude is greater than the magnitude of rate-ofstrain. So, these results validate the baffles effect on the vortex appearance inside stirred vessels.

For the baffled case baffles inhibit the vortex formation, which increase the mixing performance.

On the other hand, the deformation of the free surface in the unbaffled preneutralizer make mixing operation more difficult and less efficient. So this results validate the baffles effect on the vortex appearance inside stirred vessels. For the baffled case baffles inhibit the vortex formation, which increase the mixing performance. On the other hand, the deformation of the free surface in the unbaffled PN make mixing operation more difficult and less efficient.

Flow homogeneity and mixing quality

To investigate the mixing quality within the baffled and unbaffled preneutralizer reactor, we tracked the deformation of gas-liquid interface during the simulations. Two screen-shots of the flow contours have been token in the initial instant (t=0s), and at t=20s. The obtained CFD results proves the baffles interest in industrial STRs to ensure an efficient homogeneity in the system, especially for reacting one, since the mixing quality is crucial to control directly the production yield.

Conclusion

In that chapter, we discussed a two phase CFD model calculations of a turbulent multiphase flow carried out for the preneutralizer reactor. The effect of baffles presence on the hydrodynamic performances and the vortex appearance inside the preneutralizer was investigated. We have found that the baffles have a significant effect on the power consumption, they represent barriers to the flow which required the double of the normally consumed power. Two recirculation loops are formed for the studied cases that confirm the radial flow generated by the Pitched Blade Turbine for this specific case as mentioned in the literature [START_REF] Roustan | Mélange: Caractéristiques des mobiles d'agitation˙Techniques de l'ingénieur[END_REF]. Also, our calculation results affirm that baffles effectively break up the vortex appearance generally characterizing the unbaffled STRs. The developed preliminary CFD model is useful to understand the flow behavior, and for identifying the key parameters that affect directly the process effectiveness. In the next chapter we will focus on an advanced CFD model that includes the preneutralization chemical reaction and a structural optimization of the current ammonia spargers design.

CHAPTER 6

CFD INVESTIGATION OF HYDRODYNAMICS IN THE PRENEUTRALIZER : ANALYSIS OF THE REACTING FLOW REACTOR FOR DAP PRODUCTION

In this chapter, a three-dimensional CFD model is developed using ANSYS Fluent, allowing accurate rotational speed of the agitator on the mass transfer in the chemical reactor preneutralizer. The hydrodynamics of multiphase gas-liquid flow is numerically investigated in the preneutralizer.

Thereafter, the CFD model sensitivity to the turbulence model is investigated by replacing the RANS kε model by three models from the LES Approach. Finally, a new gas spargers type is proposed, designed and evaluated to enhance the preneutralization reactor efficiency.

Introduction

The prediction capability of CFD models has increased in importance and in accuracy, and has been a significant benefactor of advances in high-performance computing (HPC) facilities [START_REF] Ranganathan | Computational flow modeling of multiphase mechanically agitated reactors[END_REF]; Devi and Kumar (2017). In particular, CFD has involved complex coupling of physicochemical characteristics of flows, heat and mass transfer, chemical reactions, and particle-bubble-droplet dynamics. As the majority of chemical processes in real applications involves highly turbulent flows and inter-phase transport phenomena, which are intrinsically multiscale degree of complexity, CFD has been firmly established as a fundamental tool to overcoming the multiscale and multi-physics challenges. However, the predictive power of CFD in multiphase flows and transport in chemical stirred tanks is still in its infancy, as only a few studies are reported [START_REF] Cheng | Modelling and experimental investigation of micromixing of single-feed semi-batch precipitation in a liquid-liquid stirred reactor[END_REF]. Despite their widespread use, gas-liquid stirred chemical reactors have been infrequently studied. Limited studies have been carried out to understand the effect of the process operating parameters on the hydrodynamics of the reacting gas-liquid reacting flows, as rotational speed, temperature, impeller position, etc.

The preneutralizer reactor is a non-standard stirred tank reactor, equipped with a pitched blades turbine agitator that i) promotes mass and heat transfer, ii) reduces foaming, iii) and improves ammonia dispersion in the liquid phase. The latter process occurs when Ammonia gas reacts with crude Phosphoric Acid (ACP) to produce some nitrogen containing compounds Elmisaoui et al. (2020b). Two principal compounds of major importance produced in these reactions, i.e., DAP and Monoammonium Phosphate (MAP). Other species are formed but in small quantities in such a way that they do not significantly affect the total nitrogen or P 2 O 5 in the final product. The reaction between ammonia and ACP in the reactor leads to MAP or DAP depending on the ratio of N (Nitrogen) and P (Phosphorus) ratio N/P. When the ratio N/P is greater than unity, the produced MAP in the first reaction is converted to DAP in the second one [START_REF] Campbell | Ammonium phosphate slurry rheology and particle properties-the influence of Fe (III) and Al (III) impurities, solid concentration and degree of neutralization[END_REF]. The overall reactions are listed as follows

       NH 3 + H 3 PO 4 ---→ NH 4 H 2 PO 4 + Heat NH 3 + NH 4 H 2 PO 4 ---→ (NH 4 ) 2 HPO 4 + Heat (6.1a) (6.1b)
These reactions occur immediately after Ammonia gas (NH 3 ) and ACP (H 3 PO 4 ) are mixed. DAP is produced in two consecutive steps, resulting in lowering the slurry viscosity and reducing the ammonia loss. In the first step, ACP with 40 to 42 weight percent of P 2 O 5 is introduced through the and after reacting with ammonia, the precipitated product with N/P of about 1.3-1.4 is carried to the granulator drum. The second step is then completed by injection of further ammonia into the drum, which keeps the N/P of the granulated product at about 1.8-2. Recycled seed granules are fed to the granulation drum. Produced slurry in the pipe reactor is sprayed onto the seed granules, and granule growth occurs. Granules leaving the granulator are first dried and then screened to separate out the product size.

Numerical methodology

In this study, the numerical simulations have been carried out using workbench of the CFD commercial software ANSYS Fluent 2021R2.

Physical model of the preneutralizer

• 45 • 90 • 135 • 180 • 315 • 270 • 225 • 0 • 100 × 100 × 6 mm renforcing bands

Numerical model

The continuity of the CFD developed model is detailed by the following formulation, that conserve the kε turbulence model to consider the turbulent flow, and the chemical reaction, with the Eddy dissipation TCI models are also coupled. Introducing the volume fraction field φ of each phase , the governing equations for the two phases are expressed as follows

                       ∂ t φ + ∂ i u ,i φ p = 0, ∂ t φu ,i + ∂ j φu ,j u ,i = g i + ∂ j τ ,ij + F ,i ϱ p , ∂ t (φY ) + ∂ j φY u ,i = -D∂ i (φ∂ i Y ) + φR ρ , (6.2a) (6.2b) (6.2c)
where the letter t represents the time and the symbols ∂ t and ∂ i denote the partial differential operators ∂/∂ t and ∂/∂x i , respectively. Furthermore, the summation convention for repeated indices is used. Here, ϱ the density, the velocity vector is denoted by u, and the gravity by g. The sum of the volume fraction of the different phases is always unity. The stress tensor is

evaluated as τ ,ij = -∂ i p + 2φµ S ,ij -S ,kk δ ij /3 with S ,ij = ∂ j u ,i + ∂ i u ,j /2.
The last term F in Eq. (6.2b) denotes the momentum transfer between the phases, which is decomposed into drag forces F d , lift forces F l , and turbulent dispersion forces F vm and it is expressed therefore as

F ,i = F d + F l + F vm . (6.3)
In the present two-fluid calculation, the following functional forms [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF] are

retained                F d = 3 4 ϱ w d s C d ∥u w -u s ∥ (u w -u s ) , F l = ϱ w C l (u w -u s ) × (∇ × u w ) , F vm = ϱ w C vm (D w t u w -D s t u s ) , (6.4)
where D t = ∂ t + u ,j ∂ i stands for the substantive derivative corresponding to each phase , and C d , C l and C vm are the drag, lift and virtual mass coefficients, respectively. For the drag coefficient, the solid particle drag model of Schiller and Naumann (1935b) is retained, which reads

C d =          24
Re 1 + 0.15Re 0.687 , if Re ≤ 10 3 , 0.44 , if Re ≥ 10 3 , (

where Re = ϱ s ∥u r ∥d s /µ s is the Reynolds number bases on the particle diameter d s and the relative velocity u r = u w -u s . In Eq. (6.2c), Y is the local mass fraction of species in phase, D is the diffusivity, R is the reaction rate of the creation/destruction of -th species Fluent (2011).

In Eq. (6.2c), the reaction rate for the species is given by

R = ν ′ , M w,i Aρ ϵ k min R Y R ν ′ R,r M w,R (6.6)
Based on the pressure correction algorithm SIMPLE method with a value of restitution coefficient of 0.9, the flow equation (6.2) is solved by numerical method, the turbulent model which is a k-ϵ turbulence model, unsteady constants using two order implicit scheme, convection, and diffusion using two order upwind scheme and two order central scheme. The interphase mass transfer and momentum transfer was calculated implicitly. From using an implicit algorithm to discretize the species transport equation, interphase mass transfer and preneutralization chemical reaction rate source term were solved explicitly.

Mesh Grid sensitivity

ANSYS Mesher (ANSYS Inc.) was used to create the three-dimensional mesh grid of the computational domain, that we discretized into unstructured cells, with specific refined zones (cf. Fig 6 .3). A grid-dependency study was carried out to evaluate mesh suitability for the the present configuration. Three mesh featuring a number of nodes of 1066150, 1599195 and 1760629, defined as coarse, medium and fine meshes, respectively. Refinements were carried out in the shaft region, near to the nozzles, at the inlets, and at the outlet of the domain, allowing for a minimum and maximum mesh size of 30 mm and 55 mm for the coarse mesh. For the medium mesh, element sizes were obtained with a minimum and maximum mesh of 25 mm and 15 mm, whereas for the fine mesh with values of 20 mm and 5 mm, see Table . 6.2. For the fine mesh, the first mesh adjacent to the nozzles is set to a width of 0.15 mm, which results in ) that varies from 30 to 50 along the length of nozzles. it is noteworthy to mention that y + representing the structure of boundary layer is an important parameter in the meshing process since it gives significant insights into the accuracy of the simulation with low Reynolds number models to ensure y + < 5 (ideally y + < 1.0) and y + > 30 for turbulent flows with high Reynolds number [START_REF] Kern | Sensitivity analysis in 3d turbine cfd[END_REF]. The inflation of the boundary layer is defined for the different mesh grids by a transition ratio of 0.272 and layers growth rate of 1.2. The mesh is of good quality with an average skewness of 0.52, a maximum non-orthogonality of 26.5 and an average non orthogonality of 1.2. The unstructured tetrahedral mesh is set with prismatic elements to allow for a y + insensitive to the mesh refinement Manual (2009) and ensure the solution independence to the mesh grid quality. The velocity magnitude along the vertical line within the preneutralizer at x 2 = 0.35 m for different mesh sizes is portrayed in Fig. 6.4. From these results, the fine and medium mesh sizes give practically the same trend of the velocity magnitude with close values. However, the prediction of the degree of reaction for the coarse mesh is found to be far to that obtained using higher cell numbers of the meshing. This fact confirms again our meshing choice discussed above. Consequently, considering the limitations of our computational resources as well as the computational time for the solution convergence in Table 6.2, the medium mesh grid seems to be the suitable choice for predicting accurately the stirred reactor performance. The PBT turbines ensure an axial flow in standard stirred tanks, and a radial flow for specific cases depending to the geometry design modifications and the impeller rotational speed.

In the present case, the sensitivity to the rotational speed is investigated. However, this rotational speeds ensure the required conditions for the chemical species interphase exchange, with a moderate velocity magnitude in the simulation domain compared to the sparged gas inlet velocity. More than that, the recirculation loop is more intense and strong with 120 rpm rotational speed which increases the chemical species interactions and ensures the good mixing. Otherwise, based on the following table 6.3, the increase of the agitator rotational speed enhances the power of the system, normally due to the increase Reynolds number expressed as function of the rotation speed. But the power and the pumping numbers becomes very small based on operating rotational speed (68 rpm) results. However, [START_REF] Bao | Influence of the top impeller diameter on the gas dispersion in a sparged multi-impeller stirred tank[END_REF] and co-workers, more than other experimental studies highlighted the depression cavities phenomenon in the under impeller zone, when the rotational speed exceeds a specific interval depending to the impeller type. The obtained results from our CFD simulations confirm this finding. As shown in Figure 6.6, the radial illustration of the gas volume fraction in different level of the tank confirms the great effect of the rotational speed in the gas phase dispersion enhancing.

Since the dispersion is more efficient near to the impeller due to the strength of the turbulent eddies in this zone.

Based on this results, the preneutralization yield improvement can be directly controlled by the rotational speed mastering. But to counteract any adverse unexpected phenomena, the variation should be done carefully.

Effect of turbulence model on the flow hydrodynamics prediction: LES SubGrid-Scale turbulence models

For the CFD modeling studies of multiphase flows in stirred reactors, different turbulence models have been developed to predict the hydrodynamics of the flow, each one of them demonstrating its effectiveness and robustness in each particular application [START_REF] Derksen | Assessment of large eddy simulations for agitated flows[END_REF]. The age-old CFD modeling method based on the Reynolds Averaged Navier-Stokes (RANS) approach is still considered the most employed turbulence modeling approach due to its acceptable accuracy and affordable computational cost for predicting the hydrodynamics involving complex geometries, but with the development of the HPC capabilities, the employment of the hybrid Detached Eddy Simulation (DES) model, which blends the RANS approach with LES founds its applications. The DES methodology shows its increasing potential and accuracy in flow prediction more than the RANS models. Especially for its informative ability to predict the flow regimes in the moving zone characterizing the impeller rotation [START_REF] Stanly | Large-eddy simulation of a wind turbine using a filtered actuator line model[END_REF]. This sction aims at evaluating the performance of three main Large-Eddy Simulation (LES) subgrid-scale (SGS) turbulence models in the preneutralizer reactor.

Physical model of the preneutralizer

In this study, the numerical simulations have been carried out using the workbench of the CFD

                       ∂ ūi ∂x i = 0 ∂ ūi ∂t + ūi ∂ ūj ∂x i = - ∂ P ∂x i + v ∂ 2 ūi ∂x 2 j - ∂τ ij ∂x j τ ij = ūi ūj -ūi ūj (6.7a) (6.7b) (6.7c)
where τ ij is the sub-grid scale stress tensor, which reflects the effect of the unresolved scales on the resolved scales. Subgrid-scale stresses resulting from the filtering operation are unknown, and require modeling. The subgrid-scale turbulence models employ the Boussinesq hypothesis as in the RANS models, computing subgrid-scale turbulent stresses from:

τ ij = 1 3 τ kk δ ij -2µ t Sij (6.8)
Three SGS turbulence models are considered in this work:

• Smagorinsky-Lilly model

µ t = ρL 2 s | S| (6.9)
• Wall Adapting Local Eddy Viscosity (WALE) model

µ t = ρL 2 s S d ij S d ij 3/2 Sij Sij 5/2 + S d ij S d ij 5/4
(6.10)

• Algebraic Wall-Modeled LES Model

µ t = ρ min (κd wall ) 2 , (C k ∆) 2 • S • 1 -exp (-(y + /25) ) (6.11)
where d wall is the wall distance, S is the strain rate, κ = 0.41, and C k = 0.2 are constants.

To consider chemical species transport, the previous system of equations is coupled to the species transport equation, as described in [START_REF] Elmisaoui | Cfd investigation of industrial gas-liquid preneutralizer based on a bioreactor benchmark for spargers optimization[END_REF].

Numerical simulation conditions

As previously stated in many CFD simulations of stirred reactors, a variety approaches may be used to cope with the movement of the impeller blades. The Multiple Reference Frame (MRF) approach was employed in this investigation [START_REF] Elmisaoui | Cfd investigation of industrial gas-liquid preneutralizer based on a bioreactor benchmark for spargers optimization[END_REF], and the SIMPLE algorithm was employed for pressure-velocity coupling.

Computational domain

The computational grid is composed of two parts: an inner spinning cylindrical volume holding the turbine and an outer, stationary volume containing the remainder of the tank. In both portions, structured grids of non-uniformly distributed hexahedral cells are utilized. To provide a more precise description of the impeller, the grid employed in the impeller region is refined.

In the tank, there are a total of 970997 grid nodes.

It is worth noting that a grid-dependency study is carried out using ANSYS Mesher to evaluate mesh suitability for the present configuration. A three-dimensional mesh grid of the com- Three mesh featuring grids are generated and defined as coarse, medium, and fine meshes, respectively. Refinements are carried out in the shaft region, near to the nozzles, at the inlets, and at the outlet of the domain, allowing for a minimum and maximum mesh size of 25 mm and 15 mm for the medium mesh grid. The last refinement level is retained based on a mesh grid sensitivity, as the grid allows for an accurate numerical solution at acceptable time costs.

Results and analysis

To evaluate the flow hydrodynamics in a three dimensional multiphase turbulent flow, the identification of coherent structures is crucial, which allows a better understanding of the real space dynamics of turbulent flows. One of the most commonly used methods for categorizing threedimensional flow structures is the Q-criterion proposed by [START_REF] Hunt | Turbulent shear flows over low hills[END_REF]. It is defined in terms of the instantaneous velocity gradient tensor as:

Q = 1 2 (∇ • u) 2 -∇u : ∇u T = 1 2 (∇ • u) 2 + ∥Ω∥ 2 2 -∥S∥ 2 2 (6.12)
Here, S is the strain rate tensor and Ω is the rotation rate tensor.

The turbulence associated with the flow through the inlet and across the stirred reactor is visualized at two different instants t = 100 s and t = 190 s in Figs. 6.8 and 6.9. As can be noted for the three SGS-models and as the flow is introduced from the nozzles into the tank, it is clear to have higher velocities near the inlet pipes. In addition, there is a velocity difference between the nozzle stream and the rotating flow in the tank, which creates local shear zones.

As a consequence, inlet nozzles become the inherent turbulence generators. The intensity of the generated turbulence is less pronounced in the classical Smagorinsky model, while it seems to be more pronounced in the WALE and WMLES models. Furthermore, for the latter models, the instantaneous flow seems to be characterized by nearly concentric streamlines in a circular fluid motion, which is due to tangential inflow. The resulting vortex column around the central outlet structure is wrapped by the external turbulent filaments, leading to intensive momentum transport.

Figure 6.10 presents the profiles of instantaneous magnitude velocities in two locations of the domain (denoted by P1 and P2 in Fig. 6.10a).

As a considerable throughput of liquid is injected near to position P2, a high velocity peak is observed in a short time. Then, as the mixing gets enhanced, the velocity decreases gradually until t = 150 s and a considerable flow amount of Ammonia gas achieves the top of the domain.

This explains the emergence of new turbulent eddies in the P2 location. LES models were able to capture the generated eddies coming from the liquid injection and from the mixing zone.

WALE model gives higher values of the velocity magnitude than WMLES and Smagorinsky-Lilly models. However, based on the P2 probe position near to the shaft wall, the WALE is more appreciated to be used for near to the wall turbulence eddies prediction. As shown in 6.6. CFD investigation of industrial gas-liquid preneutralizer based on a bioreactor benchmark for spargers optimization

The performance of the PN depends mainly on the operating parameters and geometrical characteristics Elmisaoui et al. (2020a). The type of agitator, the design of the injectors and their positions are key parameters to improve the hydrodynamics, and thus better control the transport phenomena between the chemical species present in the system. The current geometrical conditions of the ammonia pipe injectors generate losses of a high rate of gaseous ammonia.

The hydrodynamics in the area near to the slurry outlet, and the power of the flow absorbed by the pump to the granulator, allows considerable amount of ammonia to leave the reactor just after injection. This is especially true for feeds from radial spargers near to the slurry outlet, which requires to put them in off mode in the current process and operate the PN with only six injectors, which decreases considerably the efficiency of the preneutralizer. The maturity of gas liquid chemical reactors in biotechnology, especially in oxidation bioreactors and fermenters, proved that gas sparger design governs the flow hydrodynamics in gas liquid sparged stirred reactor [START_REF] Birch | The influence of sparger design and location on gas dispersion in stirred vessels[END_REF].

It has been demonstrated that gas sparger design, dimensions and location have significant effect on the gas dispersion and the flow behavior. Many designs are used (conical sparger, jet distributor, ring sparger, tubular sparger, etc.) and the performance of each type depends to system characteristics [START_REF] Li | Review on the Effect of Gas Distributor on Flow Behavior and Reaction Performance of the Bubble/Slurry Reactors[END_REF]. Hence, the main objective of this study is to use CFD modeling to investigate the effect of changing the radial classical ammonia sparger by a new design based on ring sparger and to quantify the PN performance in terms of mixing and mass transfer.

As presented in chapter 3, the ammonia gas reaction with crude phosphoric acid (ACP)

(Liquid) to produce some nitrogen containing compounds (slurry). Two principal compounds of major importance forming in this reaction are Di Ammonium Phosphate (DAP) and Mono Ammonium Phosphate (MAP). The chemical equation showing the principal reaction is listed below:

CFD modeling strategy

The numerical model is based on the Eulerian Multi Fluid Approach, and has been carried out using the commercial CFD software ANSYS Fluent 2020 R1. Details of the model will be described in the following subsections.

Computational domain and geometry design

The PRN geometrical model is based on the industrial scale dimensions, as show in Fig. 6.11, the working volume of that reactor is 47 m 3 to ensure a throughput around 120 m3/h. It is equipped by a PBT agitator, with four impellers slanted with 45 • C from the horizontal plane. The impeller diameter is 2.5 m and is placed at a height of 0.9 m from the bottom. The impeller motion was modeled by the inner-outer model, based on the MRF approach. The ACP liquid is injected from the top, and the ammonia spargers are implemented near to the bottom.

Numerical model

The solver setting included phase coupled algorithm for pressure-velocity coupling and a finite volume method. The special discretization was set to QUICK for the volume fraction, combined with a second order scheme for the turbulence equations and the second-order backward Euler scheme for the transient terms. The under relaxation factors of the numerical scheme were set initially to 0.3 for pressure, 0.4 for momentum and 0.2 for the species volume fraction.

The convergence criterion was set to 1 × 10 4 for scalars, chemical species volume fractions, and velocities. For both cases simulations were run for a period of 50 s, and a reasonable convergence was achieved with residuals less than 10e -5 for all the equations of the system.

Domain meshing and numerical model

The ANSYS meshing tool is used to generate an unstructured grid composed by tetrahedral elements. Mesh quality has been evaluated and mesh convergence analysis has been conducted

to use the optimal grid in terms of computational time and accurate results. Figures 6.12. and 6.13 show cross-sectional and frontal views of these meshes. In this work the Euler-Euler multi-fluid model is adopted where gas and liquid phases are all treated as inter-penetrating continua, in the computational domain. The Euler-Euler (E -E)

approach is suitable for modeling such two fluid flow in a stirred reacting tank with the sparged The equation for conservation of momentum in the i th direction is defined by:

∂ ρ ϕ α ϕ ∂t + ∇ • ρ ϕ α ϕ u ϕ u ϕ = -α ϕ ∇p + ∇ • µ ϕ,ef f α ϕ ∇u ϕ + ∇u ϕ T + ρ ϕ α ϕ g + F ϕ (6.15)
Where ρ ϕ , α ϕ , and u ϕ are the density, volume fraction, and mean velocity vector of phase ϕ, respectively, where the subscript ϕ refers to the liquid (I) or gas (g) phase. The source term S m models the mass added through phase changes.

To solve the conservation equations for the chemical species in multiphase flows, within each phase q, we need to predict the local mass fraction, Y q i , through the solution of a Convection-Diffusion equation for the i th species. The chemical species conservation equation, when applied to a multiphase flow can be represented in the following form

∂ ∂t ρ q α q Y q i + ∇ • ρ q α q ⃗ v q Y q i = -∇ • α q ⃗ J q i + α q R q i + α q S q i + n p=1
ṁp i q j -ṁq j p i + R (6.16)

where R i q is the net rate of production of homogeneous species i by chemical reaction for phase q, ṁq j p i is the mass transfer source between species i and j from phase q to p, and R is the heterogeneous reaction rate. In addition, α q is the volume fraction for phase q and S i q is the rate of creation by addition from the dispersed phase plus any user-defined sources. The net rate of production of species i due to reaction r, R i,r , are previously expressed term by term in chapter 3. The energy equation is also considered to represent heat exchanges.

Flow field hydrodynamics

CFD model results show that the hydrodynamics within the preneutralizer depend directly on the spargers type. Figure 3 represents the ammonia gas distribution at the radial plane splitting the spargers from the middle, after 30s of simulation time, allowing the gas phase to attain the spargers outlets and to mix with the ACP liquid. In Figure 6.14. (a) the gas injected by the pipe sparger near to the slurry is subject to the aspiration form the pump in this boundary, which represents the same phenomenon observed in the industrial facility.

A proportion of 80% of the ammonia gas injected is well dispersed overall the PN with the new sparger design, otherwise in the zone with the principal ammonia inlet, the injected gas phase is imported to the other zone, by mixing. Also, the shape of this ammonia distributor represents an obstacle to inhomogeneous flow which impacts the hydrodynamics and transfer mechanisms further. Changes about the position of this latter are required.On other hand, the combination of a PBT and ring spargers allows a good axial dispersion of the injected ammonia Figure 6.14. (d), compares to the dispersion observed with a pipe sparger Figure 6.14. (c). The gas attains the turbulent zone very fast that for the pipe injector case, which will increase the reaction rate and decreases the required time to achieve system homogeneity.

Gas holdup distribution and mass transfer

When the ammonia gas is sparged in the PRN, it rises upwards due to the buoyancy force associated with this phase. It then diffuses or escapes in the liquid phase. We defined six probe points to follow the gas holdup within the PRN. It was observed that, for P1, P2, P3, P4 and P5 the ammonia gas is detected for both designs, and the gas dispersion decreases when rising in the vessel. At the radial plane containing P6 point (z = 4.5 m), the NH 3 phase is not detected with geometry representing pipe spargers, but it attains some specific zones with the new proposed geometry, which reflects that only some spargers are effective (Figure 6.17 (b)). Previous experimental study validated our results [START_REF] Rewatkar | gas hold-up behavior of mechanically agitated gas-liquid reactors using pitched blade downflow turbines[END_REF]. Hence, the change of the sparger type to the ring one implicates new operating conditions. For instance; the gas flow rate or geometrical characteristics (sparger diameter, location, ring diameter, etc.) are manipulated for a good hydrodynamics. This effect is more and more pronounced in the case of ring spargers when the size of injectors exceeds 0.3 m. From other point of view, the mass transfer coefficient is enhanced by increasing diameters. Figure 6.16 shows higher values with the ring sparger compared to the pipe sparger over the different diameters, but the flow pattern in the two biggest one (0.416 m -0.52 m) of this injector is quite different from those generated by the other three spargers. This is due to the stirred sparged gas liquid reactors design, where the sparging diameter and the vessel diameter

Conclusion

A 3D Euler-Euler gas-liquid CFD model is developed using ANSYS 2020 R1, allowing to evaluate the effect of sparger design on the hydrodynamics within the preneutralizer. k -ε turbulence model and MRF approach were used. The mass transfer, gas holdup and dispersion are numerically investigated. A new ring sparger is designed inspired from bioreactors, and its effect is investigated on the preneutralizer performance, especially for mass transfer and gas holdup.

The model provides good agreement with previous literature studies, and confirms our design recommendations. The new designed sparger proves its efficiency in enhancing the ammonia dispersion within the preneutralizer. The mass transfer coefficient and gas holdup sensitivity to the sparger diameter was ex animated, and the actual diameter increasing leads to enhancing these key parameters, up to the ratio d sp /D = 0.8. Further work is ongoing to improve the preneutralization manufacturing by experimenting the proposed sparger design.

In the second part of that chapter, we studied the turbulence model effects on the flow hydro-dynamics predictions and vortex detection. The liquid and gas phase flow fields in the stirred tank simulated by the LES approach are more asymmetric in the case of WMLES and WALE models. Overall, each of the three models has its own specific ability for the flow prediction.

However, the choice of an LES Turbulence model should be done carefully, for a good physical flow prediction, with a fine mesh grid, which is very costing in terms of calculations time and hardware resources compared to the RANS turbulence model, allowing to have significant results and the accurate model solution with an acceptable average of calculations time, with less requirements on the grid refinements. So, RANS turbulence model is adopted for the rest of studies. Specially, RANS k-ε model, validated the more accurate for 3D industrial scale CFD simulations.

Introduction

Rotating drums are present in many industrial processes; pharmaceutical, fertilizers, and food industries. Their areas of applications include granulation, cooling, coating, milling, reactions, drying, and mixing [START_REF] Zheng | Numerical simulation of flow characteristics in an inclining rotating kiln with continuous feeding[END_REF]. As a process, it may penalize these industries in terms of profitability, quality of finished product and respect of commercial and environmental constraints [START_REF] Mellmann | The transverse motion of solids in rotating cylinders-forms of motion and transition behavior[END_REF]. Indeed, this operation presents a serious challenge for manufacturers, because its performance depends heavily on the operating parameters and the geometrical characteristics.

CFD is used extensively to solve numerous industrial challenges due to the maturity of computational methods and advanced computing speeds and its ability to predict the results of real experiments. Different numerical approaches have been proposed to predict the transverse and axial particle motion and particles size distribution in rotary equipments. Physical models for the particles flow analysis have been developed to relate the flow behavior to the key operating parameters as rotational speed, filling level, and inclination angle.

Model development

Computational Fluid Dynamics (CFD) has proven its effectiveness as a predictive tool for complex flows modeling. Regarding multiphase CFD modeling of granular flows in rotary drums, both the Eulerian approach and the Lagrangian approach are used. Several studies Arruda et al. (2009b);[START_REF] Sm Nascimento | Solids holdup in flighted rotating drums: An experimental and simulation study[END_REF] used Discrete Element Method (DEM) as lagrangian simulation approach to study different rotary drum configurations : with and without flights. The numerical results agreed well with experimental observations, but for all this cases the simulations were done either with a reduced dimension of the drum; or the particles diameter with a specific ratio that define the number of particles in the axial direction.

The mean advantage of DEM simulations is their ability to take into account the particles shape characteristics, and allowing to the track the particles through the computational domain.

However it is limited to reduced scale equipments due to the computational time.

Alternatively, the Eulerian-Eulerian (E-E) approach is commonly adopted for CFD simulations involving several (N) phases, and these phases are assumed to form interpenetrating continua. The dynamic behavior of a multiphase system with air and fertilizer granules can be described through Multiphase Navier-Stokes (NS) equations composed by continuity and momentum equations to be solved for each phase in the Eulerian framework. Also the E-E ap-proach is less computationally expensive, and is the most widely applied approach in the study of the hydrodynamics of the multiphase flows including granular phases. Due to the continuum assumption of the particle phase,the E-E approach requires additional closure models. The kinetic theory of granular flow (KTGF) is a tool developed from the kinetic theory of gases to describe the kinematics of a granular media. It is used to describe the particle motion in fluidsolid flows due to the inherently random particle movement caused by the fluid mechanics. It (KTGF) proves its effectiveness in rotary drum systems hydrodynamic studies, when coupled to NS equations in CFD model.

Model description

To model a multiphase flow, it is essential to take into account the conservation laws and the interactions between the present phases. In this work, an Eulerian multiphase model coupled with a KTGF model is used to model the air and fertilizer granules flow inside the inclined rotary drum granulator. The model solves momentum and continuity equations for each phase.

The continuity equations for gas-solid mass conservation is given by the following equations:

             ∂ α g ρ g ∂t + ∇ • α g ⃗ v g ρ g = 0 ∂ (α s ρ s ) ∂t + ∇ • (α s ⃗ v s ρ s ) = 0 (7.1)
where, ⃗ v g and ⃗ v s are the velocities of air and fertilizer particles. The subscripts g and s stand for the gas phase (air) and the solid phase (particles) respectively. The air volume V g is defined by:

V g = I V α g dV (7.2)
Details on the air and solid momentum, particles phase stress description, collision and drag forces, and the particles characteristics from the KTGF are given in Santos et al. (2013b). The conservation of the gas and solid momentum is given by the following equations:

∂ ∂t α g ρ g ⃗ v g + ∇ • α g ρ g ⃗ v g = -α g ∇p + ∇ • -→ τ g + α g ρ g ⃗ g + α g ρ g ⃗ F g + ⃗ F lif t,g + ⃗ F vm,g + K gs ( -→ v g --→ v s ) (7.3) ∂ ∂t (α s ρ s ⃗ v s )+∇•(α s ρ s ⃗ v s ⃗ v s ) = -α s ∇p-∇p s +∇•s s +α s ρ s ⃗ g+α s ρ s ⃗ r q + ⃗ F lif t,s + ⃗ F vms + K gs ⃗ v s -⃗ v g (7.4) r0.5

Computational domain

One important step in CFD modeling is the mesh construction that represents the computational grid. A grid convergence test is done, in which the computational domain is discretised using three different grids containing a coarse, medium and fine mesh tetrahedral cells. The solid fraction distribution were determined using the different grid sizes, and the results are presented in Figure 7.2. Since the difference between numerical predictions obtained using the fine and the medium grids is negligible, the medium grid (Figure 7.1) was selected to maintain a balance between computation load and numerical accuracy. At the middle of the fluid domain by increasing the rotational speed the transition from the "slipping" to the "rolling regime", also the gravitational forces impacts more the particles flow behavior due to the slope at that level, and in the zone until the outlet of the drum a "crushing regime" characterizes the flow.

Conclusion

Solid-gas two phase particle dynamics in a rotary drum granulator was modeled and simulated using the Eulerian Granular Multiphase Model, with coupling the Euler-Euler approach to the

Introduction

In the present chapter, a CFD simulation of the multiphase flow within an industrial rotating drum granulator is carried out. The aim is to analyze the sensitivity of the rotating drum granulator to the operating parameters and to evaluate their effects on the hydrodynamics. The gas-solid two-fluid flow inside the inclined rotating drum granulator is considered and the corresponding mathematical model describing the multi-phase flow is presented in Section 8.2.

The three evaluated drag models are also presented. After presenting the physical and numerical parameters for the simulations in Section 8.3, the flow predictions are discussed in Section 8.4, focusing on the radial and axial particle distribution dependence to the operating conditions. The accuracy of the hydrodynamics prediction is analyzed and some effective parameters range of validity is proposed for the different drag models.

CFD model setup

As mentioned above, the Eulerian-Eulerian unsteady TFM is used in its general formalism, and is coupled to the kinetic theory of granular flow equations. The model conservation equations are solved using the Ansys Fluent 2020 solver.

Governing equations and fluid-solid interactions

Neglecting the effect of compressibility, the continuity equations for gas-solid mass conservation are represented by the following equations :

∇ • α g v g ϱ g = 0, ∇ • (α s v s ϱ s ) = 0, (8.1)
where, the subscripts g and s stand for the gas phase (air) and the solid phase (particles), respectively. In Eq.(8.1), α, ϱ, and v stand for the volume fraction, density and velocity, respectively, with α g + α s = 1. The conservation of the gas and solid momentum is given by the following equations (8.3) where F ext q (q stand either for g or s), F lift q , and F vm q are the external body force, the lift force and the virtual mass force, respectively. τ q are the phase stress tensors given by :

∂ ∂t α g ϱ g v g + ∇ • α g ϱ g v g v g = -α g ∇p+ ∇ • τ g + α g ϱ g g + α g ϱ g F ext g + F lift g + F vm g + K gs v s -v g , (8.2) ∂ ∂t (α s ϱ s v s ) + ∇ • (α s ϱ s v s v s ) = -α s ∇p + ∇p s + ∇ • τ s + α s ϱ s g + α s ϱ s F ext s + F lift s + F vm s + K gs v g -v s ,
τ g = α g µ g ∇v g + ∇v ⊺ g τ s = α s µ s ∇v s + ∇v ⊺ s + α s λ s -2 3 µ s ∇ • v s I (8.4)
K sg is the momentum drag exchange coefficient between the phases. In the current investigation, we assume that the system is isothermal, and only drag and gravity forces are considered.

Indeed, many research studies Santos et al. (2016a); [START_REF] Delele | Studying the solids and fluid flow behavior in rotary drums based on a multiphase CFD model[END_REF] concluded that the effects of the lift and virtual Mass forces are negligible compared to the drag force impact. We confirmed this with three simulation cases. The momentum exchange coefficient K sg between solid phase s and fluid phase g can be written in the following general form

K sg = α s ϱ s f d ϖ s ,
where f d is the drag function and ϖ s is the particulate relaxation time defined as

ϖ s = ϱ s d 2 s 18µ g ,
where d s is the diameter of particles in the solid phase and µ g is the air viscosity. The definition of the drag function f d includes a drag coefficient C d that is based on the relative Reynolds number Re r defined as

Re r = ϱ f v g -v b d s µ g . (8.5)
For the C d coefficient, different empirical and semi-empirical models exist and some will be presented in the next subsection.

Drag forces

The 8. In silico CFD investigation of the granulation hydrodynamics in rotating drum [START_REF] Wen | Chemical engineering progress symposium series[END_REF] and is given by The second one was introduced by [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF] to correct [START_REF] Wen | Chemical engineering progress symposium series[END_REF] 's model following the expression

K sg = 3C d α s α g ϱ g v s -v g 4d s α -2.
K sg =          3C d α s α g ϱ g |vs-vg| 4d s α -2.65 g for α g > 0.8, 150α s (1-α g )µ g α g d 2 s + 1.75ϱ g α s d s v s -v g for α g ⩽ 0.8. (8.8)
The third considered model is due to [START_REF] Syamlal | The derivation of a drag coefficient formula from velocity-voidage correlations[END_REF] and it can be expressed as:

K sg = 3C d α s α g ϱ g 4V 2 r,s d s Re s V r,s v s -v g , ( 8

Solid-solid interactions

To model the particle-particle interactions, the KTGF is used to define the stress term τ [START_REF] Khalilitehrani | Modeling dilute and dense granular flows in a high shear granulator[END_REF]. This theory is based on an analogy with the classical kinetic gas theory. For instance, a term θ s , representing the granular temperature, also called "pseudo-thermal", accounts for the kinetic energy associated with the particle velocity fluctuations [START_REF] Chao | Investigation of the particle-particle drag in a dense binary fluidized bed[END_REF]. The transport equation for the granular temperature is expressed as following To estimate the viscosity of the solid phase µ s , three terms are used: (8.17) where µ fr s is the frictional viscosity, µ col s is the collisional viscosity, and µ kin s is the kinetic viscosity, which are modeled using detailed correlations that can be found in [START_REF] Schaeffer | Instability in the evolution equations describing incompressible granular flow[END_REF].

3 2 ∂ ∂t (ϱ s α s θ s ) + ∇ • (ϱ s α s v s θ s ) = (-p s I + τ s ) : ∇v s + ∇ • k θ s ∇θ s -γ θ s -3K sg Φ gs , , ( 8
µ s = µ fr s + µ col s + µ kin s ,

Numerical model setup

The numerical simulations are performed for a three dimensional domain with the real industrial scale dimensions. The rotating drum granulator geometry and generated mesh are developed using ANSYS Design Modeler and ANSYS mesher, respectively. The granulator consists of a cylindrical inclined drum of diameter D = 4.5 m and length L = 9.0 m. The cylinder is inclined from the lateral plane with 3.57 • in the nominal industrial conditions. To perform the grid independence analysis, we have created three meshes with different grid refinement qualities. The meshing used tetrahedral hybrid meshes with a grid of up to 2×10 5 cells. The influence of mesh size on the rotating drum simulations will be presented in Section 8.3.2.

The second order upwind scheme is adopted to discretize the governing equations, while the QUICK scheme is employed for the volume fraction transport equation. For the pressurevelocity coupling algorithm, the pressure based Coupled SIMPLE algorithm is adopted. The transient CFD calculations are performed, and the time step size of 10 -4 s and maximum of 20 iterations per time step. We defined 10 -4 as a time step value to allow the rotary drum to rotate for more than 7 complete rotational cycles and to ensure convergence for all the rotational

Physical parameters

The primary phase is defined by air and the secondary phase is formed by the solid fertilizer particles. The particle size distribution effect is neglected, and the averaged particles diameter is set to 2 mm. The physical characteristics of the particles and air are summarized in Tab. 8.1.

The drum is initially filled with a fertilizer particles bed, and is operated in batch mode. The 1574 kg/m 3 Fertilizer particles diameter (d s ) 2 mm wall of the drum is rotating, with a rotational speed set for each simulation. It is assumed to be a non-slip wall that conveys the necessary shear traction to the granular flow as shown by Fig. 8.1. Generally, the products feed the granulator from the upper zone (inlet) and leave from the lower zone (Outlet). In this simulation, the batch mode is ensured by the definition of these boundaries as rebound periodic conditions. To investigate the fertilizer particles dynamics using the Eulerian approach, different drag models are evaluated using different operating conditions. Three different values of rotational speeds, filling degrees and inclination degrees, for the three different drag models, are evaluated as shown in Tab. 8.2. Hence, the total number of simulation cases is 81. The average CPU time is 576 core-hours per transient simulation of physical time 200 seconds. [START_REF] Wen | Chemical engineering progress symposium series[END_REF] 8.3.2. Mesh independence analysis

One of the important tasks in CFD modelling is the meshing, as there is a direct relationship between the mesh grid quality and the reliability of the simulation results. In this part of our study, three dimensional grids that describe the geometrical model and are able to handle the studied flow specifications are generated. The assessment of mesh convergence is performed with three levels of mesh refinement, as shown in Fig. 8.2.

Figure 8.2: Meshes for mesh convergence study

The computational domain is discretized using the three meshes containing 20700 (coarse), 170625 (medium) and 266045 (fine) quadrilateral cells. The two-phase simulation is considered for the grid independence/convergence analysis. For the wall layer modeling, and for the three mesh grid qualities we ensured y + in the range of 10-300 [START_REF] Chmielewski | Three-zonal wall function for k-ε turbulence models[END_REF]. The y + prediction also takes into consideration the rotating walls. An averaged value of y + < 100 in all internal domain, and y + < 11 for near the wall cells are respected in the meshing procedure.

At the initial instant (t = 0 s), the radial distribution of the two-phases is observed in Fig. 8.3.

The mesh refinement impact on capturing the interphase surface is clear. Following these results, the difference between numerical predictions obtained using the refined and the medium grids is negligible, the medium grid featuring 170625 cells is retained to keep reasonable computational costs.

Results and discussion

Mass, heat, and momentum transfer efficiency within rotary drums depends on the flow hydrodynamics. Thereby, the understanding of the flow regimes is fundamental to determine the undesirable conditions to avoid when monitoring this unit operation. For rotating drums models with a dense solid phase, the lift, virtual mass and wall friction parameters can be negligible compared to the gravity and drag model impact on the solid phase flow prediction [START_REF] Hou | Numerical and experimental study of radial segregation of bi-disperse particles in a quasi-twodimensional horizontal rotating drum[END_REF]. We confirmed this for the present computational model with three simula- prediction. In our study, the comparison results reveals that the particles velocity predictions are more accurate for [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF] and [START_REF] Wen | Chemical engineering progress symposium series[END_REF] models, than the prediction of [START_REF] Syamlal | The derivation of a drag coefficient formula from velocity-voidage correlations[END_REF] model.

However, almost identical results are obtained using the [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF] and Wen and Yu 123

Effect on the two-phase volume fraction distribution

As the drag model is a dominant element in the particles interactions prediction having the priority to be define in CFD model than other CFD subcomponents (coalescence, breakage,...), and to be combined with segregation model, for robust further CFD simulations allowing a better particles hydrodynamics prediction. In this section we aim to validate our observations on drag models effects, and their ability to accurately predict the particles phase, a projection of the drum on a transversal middle plane is considered to quantify the axial dispersion of fertilizer particles at different simulation times (see Fig. 8.6).

At the beginning of the drum rotation (t = 5 s), the interface between air and solid bed conserves its initial flat pattern. Afterwards, at t = 10 s, a portion of the particles reaches the top of the drum, and a fraction of particles adheres to the wall of the drum. On the other hand, due to the drum's inclination, most particles are stacked in the lower zone and do not reach the top of the drum. The gravitational forces acting on particles in the axial direction dominate the velocity axial component and most particles accumulate near the bottom of the drum.

At t = 20 s, most particles from the active zone attain the higher part of the drum and an Sshape surface appears on the cross section plane with significant curvature announcing a transition from the rolling to the cascading regime. Considering the effect of rotation, the bed is more expanded, which increases the action of inertial forces on the fertilizer grains. However, granular layers with lower density and higher surface velocity favor regime transition. From

Flow sensitivity to operating conditions

Rotational speed effect

In this section, the effects of rotational speeds on particle dynamics are investigated. Rotational speed is a key parameter for controlling the particle dynamics inside the granulator . The variations of the particles' average velocity at the x-y middle plane of the drum is portrayed in Fig. 8.7, for each drag model, for three rotational speeds of 8, 16, and 24 rpm, and for three filling levels of 20%, 30%, and 40% . It can be observed that the particle's velocity decreases with increasing bed depth. On the other hand, it considerably increases when increasing the rotational speed of the rotating drum. This is due to the transfer of the energy from the drum to the granular bed. The overall particles distribution within the granular bed can be explained by the fact that the thickness of the active layer is more considerable when the rotational speed is higher and the bed depth is lower, so that fertilizer particles attain the highest mean velocity with a δ of 0.4 m (filling rate of 20%), and a maximal rotational speed of 24 rpm. Increasing simultaneously the rotational speed and the depth of the granular bed leads to a decrease of the velocity of fertilizer particles. This is due to a higher thickness of the passive layer in the granular bed, which reduces the fraction of particles in the active layer, and slows particle movement within the drum. The present findings agree well with the study performed by [START_REF] Wu | Numerical study of size-driven segregation of binary particles in a rotary drum with lower filling level[END_REF].

Moreover, at the different rotational speeds and fill levels, the predictions given by [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF] and [START_REF] Wen | Chemical engineering progress symposium series[END_REF] models are quite similar. [START_REF] Syamlal | The derivation of a drag coefficient formula from velocity-voidage correlations[END_REF] model gives a prediction having a similar trend at low bed depths, but deviates when increasing the bed depth. We can conclude from these section and subsection §4.1, that although we validate the observations from previous works [START_REF] Upadhyay | An assessment of drag models in Eulerian-Eulerian CFD simulation of gas-solid flow hydrodynamics in circulating fluidized bed riser[END_REF], that recommend this drag model in non dense gas-solid systems such as fluidized beds, we point out the limitation of this models for dense solid systems such rotating drums. [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF] and [START_REF] Wen | Chemical engineering progress symposium series[END_REF] drag models are more accurate for gas-solid (dense) flows in rotating drums.

Filling degree effect

The rotating granulator is set to different rotational speeds in order to analyze the different regimes of solid motion. The results in Fig. 8.9 show the transition between the typical transverse bed motions (rolling, cascading, and cataracting) during the rotation of the rotating drum.

The configurations in Fig. 8.9 represent the hydrodynamics in a radial section in the middle of the rotary granulator for three filling rates : 20%, 30%, and 40% of the solid granular phase. We present the results for a rotational speed of 8 rpm which is the nominal speed in the industrial important deformations in its active layer for filling rates of 20 and 30%, but less deformations for the level 40%. In this latter case, the bed keeps a sliding flow regime during the simulation, hindering therefore good mixing in the granular bed and inhibiting the exchanges of mass, heat, and momentum between the granules. The case of a 20% and 30% filling rate allows the granular phase to easily follow the rotational movement of the granulator and the active layer is thicker. Indeed, a total adhesion to the walls of the domain is observed around the wall of the granulator, until reaching the height threshold causing the avalanche phenomenon, and the granular bed is put back at the bottom of the domain by gravity effect.

This effect is less accentuated when using a value of solid volume fraction of 20%, where resistance to gravity is clearer, and the granular bed is more compact, whereas the critical flow regime of cataracting is reached. This is probably due to the fact that stronger frictional forces are formed between the particles, which leads to an increase in the cohesive force that can cause a transition from the rolling regime to the cataracting regime or to the centrifuging regime with higher rotational speeds. This is not favorable in the case of a reactive flow, for which certain exchange constraints between the phases must be respected. The present results define a range of variation of filling rate to be respected with a minimum of 20% and a rate not to be reached of 40%. We note that following the present results, the rolling regime is the most adequate to achieve good hydrodynamics within the granulator.

Inclination degree effect

The tailored numerical model was employed to investigate the effect of the inclination degree on the flow hydrodynamics along the rotating drum. The impact of inclination on the mean axial velocity variation along the drum is less significant with the increase of rotational speed. Our results are consistent with Widhate et al.

(2020) to validate the impact of inclination degree on the particles dynamics within the rotating drum. Moreover, the variation of axial velocity of particles along the rotating drum explains the regimes transitions phenomena previously mentioned, especially; rolling-cascading regimes by increasing the rotational speed, with a 30% fill level. Similar observations were highlighted by the experimental study of [START_REF] Li | Axial transport and residence time of msw in rotary kilns: Part I. experimental[END_REF] and proved in experimental and simulation studies by [START_REF] Rogers | A Monte Carlo method for simulating dispersion and transport through horizontal rotating cylinders[END_REF]. Part IV

Conclusions and perspectives

to study carefully the complex flow hydrodynamics under the purpose of evaluating the effects of the structural and operational characteristics of this unit on its performance. Its worth noting that, to our knowledge; this is the first research work investigating a 3D CFD modeling on industrial scale dimensions in phosphate fertilizer industry.

We proved that CFD approach is a reliable tool to better understand multiphase flows without being under the constraint of running multiple pilot or industrial experimental tests to assess the impact of operating parameters. And with no remaining to invest huge budgets in the conception and implementation on new design to evaluate their impacts on the preneutralizer performance.

We investigated the multiphase flow considering two continuous phases :Phosphoric Acid & Ammonia, and we established a sensitivity analysis of the power consumption, the pumping capacity and the solid dispersion by varying the impeller rotational speed and reactants throughputs, more than the structural design impacts on the flow hydrodynamics. We found that increasing the impeller speed improves the pumping capacity and mixing well in the preneutralizer but it engenders an excessive consumption of energy and from a specific range of rotational speed, a depression zone have place near to the wall. So these findings were crossed to conclude the importance of the determination of the appropriate impeller speed for each manner in order to keep a great performance of the preneutralization without consuming energy excessively and avoid any depression phenomenon against the efficient hydrodynamics.

Then, the sensitivity to the operating parameters and structural details of the preneutralizer were evaluated: the two major structural modified elements are: baffles and the design of ammonia gas spargers. These, for each case the speed of the hydrodynamics performance indicators were determined, thus a better assessment were established concerning the hydrodynamic behavior of the flow within the preneutralizer reactor, and for every rotational speed, an optimal impeller speed was defined. As a result, the mixture homogenizes optimally without consuming additional energy. The quantitative study was endorsed by a qualitative study in order to identify the vortex appearance depending on the impeller speed.

The studies held in this thesis work are important for industrial evaluation of the agitation processes. We took profit of the CFD tool to investigate the hydrodynamic of multiphase flows within mechanically agitated vessels, and the CFD model efficiency link to the turbulence model namics. The results confirmed the CFD model ability to assess multiphase flow for industrial scale rotary drum systems. In light of the present work, new operating conditions combinations are proposed for the industrial plant (for example : 16 rpm, 30%, 3.57 • ).

Contributions of the thesis

The scientific findings and novel aspects of the current study are as follows :

1. The introduction of a new method for Simulating an industrial level unit operations with multiphase flows (gas-liquid and gas-solid); in three dimensions.

2. The systematic investigation of the effects of operating conditions and structural design modifications on the extent of turbulent homogeneities in the stirred chemical reactor equipped with an baffles and pitch blade turbine;

3. The proposed new preneutralizer structure in this work may be employed with limited changes for further evaluation (gas sparger's number, the specific diameter, and distance between the spargers and the agitator), and conduct to a new efficient preneutralization chemical reactor invention based on the existing process modification;

4. The development of a complete workbench of a CFD model to serve as a tool for the study of the rotating drum granulator sensitivity to the main operating conditions: Filling level, inclination degree, and rotational speed; also an advanced study is done to determine the accurate drag model as an essential subcomponent to the CFD model.

General aspects of the CFD model and perspectives

The turbulent flow was solved with a two-equation RANS approach. The realizable k-ε gave overall accurate results in terms of global mixing time, power draw, and liquid velocity, in both single and multiphase flows. Moreover, the use of periodic boundary conditions together with the MRF model, did not affect the accuracy of the model while keeping the computational cost bearable.

Simulations of gas-liquid systems were accomplished by using the Eulerian model; the interaction between the phases was described in terms of drag coefficient and was both numerically stable and sufficiently precise.

An accurate estimation of the average gas holdup, k La , and liquid velocity was obtained from the two-phases simulations. However, in some simulations, the spatial distribution of gas fraction appeared to be more heterogeneous than what has been reported in the literature.

The proposed methodology for the investigation of the multiphase turbulent chemical flow in STRs was used for the original implemented turbine (PBT). It would be interesting to try this strategy to study a multistage stirred tank, for different type of turbines. Its a promising modification that deal with the yield enhancing of complex flows in STRs.

The results of that research work constitute a data base generator to a 0D PBM internal implemented model to model the granulation process loop with a powerful industrial and academic simulator to be used for the fertilizer production stuff training.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Model development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.2 Numerical simulation conditions . . . . . . . . . . . . . . . . . . . . . . . . 7.2.3 Computational domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Results and analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.1 Flow behavior and regimes transition for air-granules simulation . . . . . . 7.3.2 Transverse flow behavior with particles size distribution (PSD) . . . . . . . 7.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 In silico CFD investigation of the granulation hydrodynamics in rotating drum 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 CFD model setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2.1 Governing equations and fluid-solid interactions . . . . . . . . . . . . . . . 8.2.2 Drag forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2.3 Solid-solid interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Numerical model setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3.1 Physical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3.2 Mesh independence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv CONTENTS xv 8.4.1 Drag model effect on hydrodynamics predictions . . . . . . . . . . . . . . . Effect on particles' velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . Effect on the two-phase volume fraction distribution . . . . . . . . . . . . . 8.4.2 Flow sensitivity to operating conditions . . . . . . . . . . . . . . . . . . . . Rotational speed effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Filling degree effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inclination degree effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Contributions of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3 General aspects of the CFD model and perspectives . . . . . . . . . . . . . . . . . . turbulence models . . . . . . . . . . . . . . . . . . . . . . . . . . . . of Di Ammonium Phosphate fertilizer affected by caking[START_REF] Becker | Phosphates and phosphoric acid: raw materials, technology, and economics of the wet process[END_REF] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Methodological strategy for CFD modeling in process engineering . . . . . . . . . 2.1 Flow Chart of DAP manufacturing process with raw material transformation. . . . 2.2 A schematic diagram of the DAP manufacturing process studied in the thesis (Technologies JESA Office, 2020). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Preneutralizer chemical reactor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Rotary drum Ammoniator-granulator . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Dryer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Classification unit: screens and crushers . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Conditioning unit: Fluidized bed cooler and rotary coater . . . . . . . . . . . . . . 3.1 Multiple reference frames representation in a baffled tank. (Cullen, 2009) . . . . . 3.2 Fixed and rotational reference frames . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Mesh types to use in computational grids development (Marshall and Bakker, 2003) 3.4 Fixed control volume for conservation laws balances . . . . . . . . . . . . . . . . . 3.5 The SIMPLE Algorithm calculation steps (Hussin et al., 2018) . . . . . . . . . . . . 3.6 Levels of turbulence models application (Marshall and Bakker, 2003) . . . . . . . 3.7 Representation diagram of two-film theory of Lewis-Whitman . . . . . . . . . . . 3.8 A schematic presentation of gas sparging in liquid phase within STR Doran (2013).

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 A schematic presentation of a horizontal and radial sections of the baffled and unbaffled standard STRs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Wet granulation technologies. Inspired from Saleh and Gugon (2009) . . . . . . . 4.2 Rotating drum granulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 The first granulation mechanism: wetting and nucleation . . . . . . . . . . . . . . 4.4 A scheme of particles growth mechanism :particle layering and binary particle agglomeration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Particle attrition and breakage illustration . . . . . . . . . . . . . . . . . . . . . . . 4.6 Regimes classification of flow configurations (Mellmann, 2001) . . . . . . . . . . . 5.1 Results of the design of the preneutralizer of the 107A production line . . . . . . 5.2 Mesh of unbaffled (left) and baffled (right) simulation domains . . . . . . . . . . . 5.3 Streamlines inside the unbaffled and baffled PN . . . . . . . . . . . . . . . . . . . . . . 5.4 Vortex inside the unbaffled and baffled PN . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Mixing development in the baffled PN (t=0s and t=20s, respectively) . . . . . . . . . . . 5.6 Mixing development in the unbaffled PN (t=0s and t=20s, respectively) . . . . . . . . . . 6.1 CAD Geometry design of the preneutralizer . . . . . . . . . . . . . . . . . . . . . . 6.2 Geometrical characteristics of the preneutralizer . . . . . . . . . . . . . . . . . . . 6.3 Mesh of the domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Velocity profiles along a vertical line (at x 2 = 0.35m) within the preneutralizer for different meshes quality (left). Area occupation ratio for different NH 3 gas distributions in the radial cross-section x 3 = 2.5m (right) with different time steps 6.5 Streamlines results projection in the middle of the domain for 68 rpm (left) and 120rpm (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Mass transfer profile in different levels (Z=3.7m; Z=2.7m; Z=1.4m) for 68 rpm and 120 rpm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Mesh grid of the simulation domain, with zoom in shaft axis and gas spargers . . 6.8 Turbulence structures for (left to right) Smagorinsky-Lilly model, WALE model and WMLES model, respectively at t = 100 s . . . . . . . . . . . . . . . . . . . . . . 6.9 Turbulence structures for (left to right) Smagorinsky-Lilly model, WALE model and WMLES model, respectively at t = 190 s . . . . . . . . . . . . . . . . . . . . . . 6.10 (a) Axial projection plan in the middle of the domain, (b) & (c) Temporal evolution of velocity magnitude profiles in probes P1 and P2, respectively. . . . . . . . . . . . . . . . 6.11 The PRN geometry with gas spargers modification, (a) original, (b) new sparger . xviii List of Figures xix 6.12 Radial cross section at the level of middle spargers plane, with zoom in the new sparger (1076818 cells) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Frontal view of the mesh domain at the middle of the PRN (1162588 cells) . . . . 6.14 Ammonia gas distribution in the radial plane splitting spargers (X=0.3m, t=40s), respectively for (a) pipe and (b) ring spargers, with a zoom on the axial distribution of ammonia in injection zone (c) pipe, (d) ring spargers . . . . . . . . . . . . . 6.15 Comparison of mass transfer coefficient obtained for pipe and ring spargers for different injectors diameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.16 Comparison of gas holdup for pipe and ring spargers for different injectors diameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.17 (a) frontal view of the PN with probes, (b) upper plane (P6) . . . . . . . . . . . . . 7.1 Mesh of the computational domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Cross sections from grid independence test . . . . . . . . . . . . . . . . . . . . . . 7.3 The regimes transitions for different rotational speeds . . . . . . . . . . . . . . . . 7.4 Solid volume fraction at the initial state . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Particles distribution along the rotary drum (particles with d 2 ) . . . . . . . . . . . . . . 8.1 Respectively, the CAD geometry and Mesh of the simulation domain . . . . . . . . 8.2 Meshes for mesh convergence study . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Radial cross section at initial time for the different meshes. . . . . . . . . . . . . . 8.4 Velocity magnitude of particles (m/s) as function of distance (a) from the wall and (b) zoom in close to the wall for the case (8 rpm, 30% filling degree, and 3.57 • inclination degree). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Contents 1 . 1

 11 Industrial context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Study objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Scientific challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Manuscript structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1. Industrial context Morocco's phosphates and derivative products exports almost doubled as of March 2022, compared to the same period in 2021 when the number reached 13.4 billion MAD. The word bank

  multiphase flow, by coupling Navier-Stokes Equations (NSE) to Turbulence Chemistry Interactions suitable models; 2. Use the previous CFD model to evaluate the sensitivity of the preneutralization performance to the operating conditions and geometrical structure of the domain design; 3. To develop numerical CFD simulations of the complex multiphase flows within a DAP rotating drum granulator (RDG) to understand the hydrodynamics within it. These CFD simulations will allow to produce the reduced models useful for process control and optimization; and to understand the particles dynamics effect on the hydrodynamics of the fluid flow within RDG, and classify the particles flow regime, and for an evaluation of the drag model subcomponent effects on the CFD model robustness in the granular multiphase flow dynamics prediction;

  Chapter 4: A literature review on the theoretical aspects of granulation, the granulation technologies, and flow regimes classification theory in rotating drum granulators. The hydrodynamics characteristics, and mass, and heat transfer mechanisms are presented in this part; • The second part presents the developments and results which have been obtained in modeling and simulation of the preneutralization; -Chapter 5: A work that focuses on the baffles effects on the hydrodynamics within a simplified preneutralizer model using CFD approach; -Chapter 6: CFD investigation of hydrodynamics in an industrial preneutralizer, in which we analyzed the reacting flow hydrodynamics and the flow sensitivity to the agitator rotational speed; and finally, a CFD based preneutralizer geometric optimization inspired by an industrial bioreactor benchmark. • The third part contains the CFD analysis of the granulation modeling and simulation, with the following chapters, -Chapter 7: A study of the fertilizer particles distribution impact on the multiphase flow behavior within the inclined rotating drum granulator; -Chapter 8: The main goal of this chapter is to assess the sensitivity of the multiphase CFD model with a granular phase to the drag model. Syamlal & O'Brien, Gidaspow and Wen & Yu models are tested. We evaluated the simulations sensitivity to the rotational speed, inclination degree, and fill level. The obtained results confirmed the CFD model ability to assess multiphase flow for industrial scale rotary drum systems.

Figure

  Figure 1.2: Methodological strategy for CFD modeling in process engineering

Figure 2 . 1 :

 21 Figure 2.1: Flow Chart of DAP manufacturing process with raw material transformation.

  Fig 2.2).

Figure 2

 2 Figure 2.2: A schematic diagram of the DAP manufacturing process studied in the thesis (Technologies JESA Office, 2020).

Figure 2 . 6 :

 26 Figure 2.6: Classification unit: screens and crushers

Figure 2 . 7 :

 27 Figure 2.7: Conditioning unit: Fluidized bed cooler and rotary coater

Figure 3 . 1 :

 31 Figure 3.1: Multiple reference frames representation in a baffled tank.[START_REF] Cullen | Food mixing: Principles and applications[END_REF]

Figure 3

 3 Figure 3.2. The rotating frame has the angular velocity ω. The position vector r from the origin of the rotating frame locates any arbitrary point in the fluid domain. The fluid velocities are

Figure 3 .

 3 Figure 3.4 represents a finite control volume defined in a flow field with streamlines. It is fixed in space and bounded by a closed surface ∂V . A surface element dS is considered as well

Figure 3 . 5 :

 35 Figure 3.5: The SIMPLE Algorithm calculation steps[START_REF] Saeed Hussin | Computational analysis of two-dimensional wing aeroelastic flutter using Navier-Stokes model[END_REF]

  Figure 3.6 presents the different turbulence models presentation.

Figure 3 . 6 :

 36 Figure 3.6: Levels of turbulence models application[START_REF] Marden | Computational fluid mixing. Handbook of industrial mixing: science and practice[END_REF]

  the most common types of turbulent models. Models like the k-ϵ and k-ω model have become industry standard models and are commonly used for most types of engineering problems. Two equation turbulent models are also very much still an active area of research and new refined two equation models are still being developed. Most often one of the transported variables is the turbulent kinetic energy, k , The second transported variables varies depending on what type of two equation models it is. Common choices are the turbulent dissipation, ε, or the Specific turbulent dissipation rate, ω. The second variable can be used to determine the scale of turbulence ( length scale or time scale), whereas the first variable, k, determines the energy in the turbulence.

  Further equations are needed to solve the system. k α k = 1. (3.22) 3. Theoretical aspects of stirred reactors: hydrodynamics, mass transfer and numerical modeling k f k = 0. (3.23) In general, Eq. (3.20) is solved for the dispersed phases, whereas the volume fraction of the continuous phase is obtained from Eq. (3.22). For the two-equations turbulence models, ANSYS Fluent allows to choose among three different class of turbulence models for multiphase flows: the mixture model, the dispersed model and the per-phase model, in which turbulent quantities are solved for each phase. In this research study, the dispersed formulation was used, because of the high concentration of the dispersed phase and the density difference between gas and liquid phase. With the dispersed method, turbulence quantities for the continuous phase are predicted from transport equations (Eqs. (3.15)-(3.18)), properly modified for taking into account the phase fraction, while turbulence quantities for the dispersed phase are obtained (Fluent, 2013).

3. 4 .

 4 Mass transfer in stirred sparged reacting flows 3.4.1. Overview on mass transfer in gas-liquid flows

Figure 3 . 7 :

 37 Figure 3.7: Representation diagram of two-film theory of Lewis-Whitman

Figure 3

 3 Figure 3.8: A schematic presentation of gas sparging in liquid phase within STR Doran (2013).

Figure 3 . 9 :

 39 Figure 3.9: Evolution at increasing gas flow rate F g of gas cavity structure behind the blades of a Rushton turbine. Original imagine from Nienow (1998).

Figure 3 .

 3 Figure 3.10: Interaction between turbulent eddies and fluid particles: a) interaction with several eddies, b) interaction with outer edge of eddy, c) transport with large eddy d) collision with paired eddies[START_REF] Andersson | Modeling the breakup of fluid particles in turbulent flows[END_REF]

Figure 3 .

 3 Figure 3.11: Turbulence Chemistry Interactions models matrix, from(Manual, 2009)

Figure 3 . 13 :

 313 Figure 3.13: Relative power demand for a PBT inclined with 45 • from[START_REF] Mmcg Warmoeskerken | Flow maps for rushton turbines[END_REF]

Figure 3 .

 3 Figure 3.14: A schematic diagram Power number as a function of Reynolds number for the main turbine designs, from[START_REF] Robert | Encyclopedia of physical science and technology[END_REF]

  .44) The pumping rate is calculated according to various vertical planes in the mixing tank. From the different pumping values withdrawn, an overall average value of the flow is obtained. The Equation (3.44) represents the surface swept by the scraping stirrer, H the height of the fluid and v + z the positive axial velocity of the fluid. Although this equation gives only the axial flow, it is also possible to adapt it to obtain the radial and tangential flow rates generated in the stirred chemical reactor. The normalization of the pumping flow rate by n/d 3 allows to quantify the pumping flow number investigating the impeller capacity to pump the fluid within the whole domain of the vessel. Its usually defined as:

  3.49)Fluid circulation is the curvilinear integral of velocity on a closed curve[START_REF] Ae Leanhardt | Coreless vortex formation in a spinor bose-einstein condensate[END_REF] to Stokes' theorem, to the surface integral of the rotational velocities on the free surface Γ = S (∇xv) dS (3.51)

Figure 3 .

 3 Figure 3.16: A schematic presentation of a horizontal and radial sections of the baffled and unbaffled standard STRs.

AFigure 4 . 1 :

 41 Figure 4.1: Wet granulation technologies. Inspired from[START_REF] Saleh | Mise en oeuvre des poudres: Techniques de granulation humide et liants[END_REF]

Figure 4 . 2 :

 42 Figure 4.2: Rotating drum granulator

Figure 4 . 3 :

 43 Figure 4.3: The first granulation mechanism: wetting and nucleation

Figure 4 . 4 :

 44 Figure 4.4: A scheme of particles growth mechanism :particle layering and binary particle agglomeration

Figure 4 . 5 :

 45 Figure 4.5: Particle attrition and breakage illustration

approach

  [START_REF] Nascimento | Eulerian study of flights discharge in a rotating drum[END_REF], Euler-Lagrange approach, Discrete Element Method (DEM) Santos et al. (2016a); Kasper et al. (2021), and Population Balance Model (PBM) are the frequent methods used for dense granular flows studies. Coupled CFD-DEM, CFD-PBM, DEM-PBM Baba et al. (2021); Zhao and Shan (2013), and CFD-DEM-PBM models are developed for a better prediction and understanding of the hydrodynamics and phases' interactions Hayashi et al. (2020).

  multiphase CFD simulation. Several research works[START_REF] Delele | Studying the solids and fluid flow behavior in rotary drums based on a multiphase CFD model[END_REF] have shown that the drag model is the critical closure model to accurately predict the hydrodynamics in multiphase gas-solid flows. However, none of these studies analyzed the flow behavior in a large scale, rotating and three-dimensional model, such as the granulation drum. Moreover, the range of validity for each model should be assessed. The accurate numerical prediction of the hydrodynamics can then be used for the optimization of the process key operating parameters on the light of the regime transition theory[START_REF] Arntz | Granular mixing and segregation in a horizontal rotating drum: A simulation study on the impact of rotational speed and fill level[END_REF];[START_REF] Pathmathas | Granular flow modelling of rotating drum flows using positron emission particle tracking[END_REF].

  The principal concept of PBM can be written by the rate of change of the density function inclass

  .22) For the C d coefficient, different empirical and semi-empirical models exist and some will be presented in the next subsection. The principal closure model for the momentum exchange equation is the drag force. The literature reports numerous Drag models employed in CFD simulations of multiphase systems containing a granular phase. Here we aim to investigate the influence of the drag model on the computation of the granular flow in a rotating drum granulator. Several research works on drag models in fluid-solid systems, give a recommendation of drag models depending on the density of the granular phase Wang et al. (2020); Du et al. (2006); Upadhyay et al. (2020).Schiller and Naumann (1935a) and[START_REF] Helland | Drag reduction and cluster formation in a circulating fluidised bed[END_REF] models are recommended for thin and diluted systems. However, for multiphase systems with dense granular phase, three models are preferred and are retained for further investigation in the present study.

Figure 4 . 6 :

 46 Figure 4.6: Regimes classification of flow configurations[START_REF] Mellmann | The transverse motion of solids in rotating cylinders-forms of motion and transition behavior[END_REF]

Figure 5 . 1 :

 51 Figure 5.1: Results of the design of the preneutralizer of the 107A production line

Figure 5 . 2 :

 52 Figure 5.2: Mesh of unbaffled (left) and baffled (right) simulation domains

Figure 5 . 3 :

 53 Figure 5.3: Streamlines inside the unbaffled and baffled PN

Figure 5 . 4 :

 54 Figure 5.4: Vortex inside the unbaffled and baffled PN

Figure 5 . 5 :

 55 Figure 5.5: Mixing development in the baffledPN (t=0s and t=20s, respectively)

Figure 5 . 6 :

 56 Figure 5.6: Mixing development in the unbaffledPN (t=0s and t=20s, respectively)

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 6.2 Numerical methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.2.1 Physical model of the preneutralizer . . . . . . . . . . . . . . . . . . . . . 85 6.2.2 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.3 Mesh Grid sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.4.1 Flow field in multiphase preneutralizer . . . . . . . . . . . . . . . . . . . 90 6.5 Effect of turbulence model on the flow hydrodynamics prediction: LES SubGrid-Scale turbulence models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.5.1 Physical model of the preneutralizer . . . . . . . . . . . . . . . . . . . . . 92 6.5.2 Mathematical modeling and numerical schemes . . . . . . . . . . . . . . 92 6.5.3 Results and analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.6 CFD investigation of industrial gas-liquid preneutralizer based on a bioreactor benchmark for spargers optimization . . . . . . . . . . . . . . . . . . . . . 95 6.6.1 CFD modeling strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.6.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

  Figure 6.1: CAD Geometry design of the preneutralizer

Figure 6 . 2 :

 62 Figure 6.2: Geometrical characteristics of the preneutralizer

  Figure 6.3: Mesh of the domain

Figure 6 . 4 :

 64 Figure 6.4: Velocity profiles along a vertical line (at x 2 = 0.35m) within the preneutralizer for different meshes quality (left). Area occupation ratio for different NH 3 gas distributions in the radial crosssection x 3 = 2.5m (right) with different time steps

Figure 6 . 5 :

 65 Figure 6.5: Streamlines results projection in the middle of the domain for 68 rpm (left) and 120rpm (right)

  commercial software ANSYS Fluent 2021-R2. The physical model tailored in this work is based on the actual size of a preneutralization chemical reactor (PN) with a working volume of 47 m 3 to ensure a throughput of around 120 m 3 /hr. The ACP is injected from the top, and eight uniform Ammonia inlet nozzles are arranged in the bottom. The rotational agitator is equipped with a Pitch Blade Turbine (PBT), consisting of four impellers slanted with an angle of 45 • from the horizontal plane. Detailed characteristics of the CAD file (cf. Fig. 7.1) are based on the original preneutralizer geometry as described in details in the previous section.

Figure 6 . 7 :

 67 Figure 6.7: Mesh grid of the simulation domain, with zoom in shaft axis and gas spargers

  6. CFD investigation of hydrodynamics in the preneutralizer : Analysis of the reacting flow reactor for DAP production putational domain is discretized into unstructured cells with specific refined zones (cf. Fig.7.1).

Figure 6 . 8 :

 68 Figure 6.8: Turbulence structures for (left to right) Smagorinsky-Lilly model, WALE model and WM-LES model, respectively at t = 100 s

  Figure 6.10: (a) Axial projection plan in the middle of the domain, (b) & (c) Temporal evolution of velocity magnitude profiles in probes P1 and P2, respectively.

Figure 6 .

 6 Figure 6.11: The PRN geometry with gas spargers modification, (a) original, (b) new sparger

  Figure 6.11 (a) and (b) represents the used geometries in this work, with the original and the new proposed spargers. In the design of the new sparger, the diameter of the ring is determined based on the distance between the original radial injectors outlets and the tank wall. This ensures the same flow throughput by each injector similarly to the original feeds. The system boundary conditions and physical settings are based on the operating conditions at the industrial scale.This reaction occurs immediately when NH 3 and H 3 PO 4 are mixed. Depending to the mole ratio between these reactants, the intermediate production of MAP happens, then the excess ammonia converts MAP ((NH 4 ) H 2 PO 4 ) to DAP ((NH 4 ) 2 HPO 4 )

Figure 6 . 12 :

 612 Figure 6.12: Radial cross section at the level of middle spargers plane, with zoom in the new sparger (1076818 cells)

Figure 6 . 13 :

 613 Figure 6.13: Frontal view of the mesh domain at the middle of the PRN (1162588 cells)

  Figure 6.14: Ammonia gas distribution in the radial plane splitting spargers (X=0.3m, t=40s), respectively for (a) pipe and (b) ring spargers, with a zoom on the axial distribution of ammonia in injection zone (c) pipe, (d) ring spargers

Figure 6 . 15 :

 615 Figure 6.15: Comparison of mass transfer coefficient obtained for pipe and ring spargers for different injectors diameters

  Figure 6.16: Comparison of gas holdup for pipe and ring spargers for different injectors diameters

Figure 6 .

 6 Figure 6.17: (a) frontal view of the PN with probes, (b) upper plane (P6)

Figure 7 . 1 :

 71 Figure 7.1: Mesh of the computational domain

Figure 7 . 2 :Figure 7 . 4 :

 7274 Figure 7.2: Cross sections from grid independence test

Figure 7 . 5 :

 75 Figure 7.5: Particles distribution along the rotary drum (particles with d 2 )

  principal closure model for the momentum exchange equation is the drag force. The literature reports numerous drag models employed in CFD simulations of multiphase systems containing a granular phase. Here we aim to investigate the influence of the drag model on the computation of the granular flow in a rotating drum granulator. Several research works on drag models in fluid-solid systems, give a recommendation of drag models depending on the density of the granular phase[START_REF] Wang | Mesoscale drag modeling: A critical review[END_REF];[START_REF] Du | Computational fluid dynamics (CFD) modeling of spouted bed: Assessment of drag coefficient correlations[END_REF];[START_REF] Upadhyay | An assessment of drag models in Eulerian-Eulerian CFD simulation of gas-solid flow hydrodynamics in circulating fluidized bed riser[END_REF].Schiller and Naumann (1935a) and[START_REF] Helland | Drag reduction and cluster formation in a circulating fluidised bed[END_REF] models are recommended for thin and diluted systems. However, for multiphase systems with dense granular phase, three models are preferred and are retained for further investigation in the present study. The first one was developed by

Figure 8 . 1 :

 81 Figure 8.1: Respectively, the CAD geometry and Mesh of the simulation domain

Fig. 8 .

 8 Fig. 8.4 presents the velocity magnitude over the perpendicular line to the particles bed in the simulation domain (xy plane, y = 2.5 m, z = 4 m).

Figure 8 . 3 :Figure 8 . 4 :

 8384 Figure 8.3: Radial cross section at initial time for the different meshes.

  tion cases. However, the global CFD model is sensitive to the drag modelSantos et al. (2013a);[START_REF] Delele | Studying the solids and fluid flow behavior in rotary drums based on a multiphase CFD model[END_REF];[START_REF] Adnan | Validation and sensitivity analysis of an eulerian-eulerian two-fluid model (TFM) for 3D simulations of a tapered fluidized bed[END_REF]. In the following, an analysis of the flow hydrodynamics and its sensitivity to drag models, rotational speed and filling level is performed to deepen our understanding of the dynamic behavior of granular flow in rotating drums.8.4.1. Drag model effect on hydrodynamics predictions Effect on particles' velocityThe time averaged results obtained on the line located at a radial cross section in (x-y plane, with y = 2 m and z = 4.5 m) of the drum is shown in Fig.8.5. The operating conditions are 30% of fill rate (that is equivalent to 0.675 m as granular bed depth), 3.57 • inclination degree from the reference plan, and a rotational speed of 8 rpm. For the three drag models, the fertilizer particles velocity profile has its highest value of 2.4 m/s at the bed surface, and then decreases with decreasing radial position to reach its lowest value at the active-passive zones interface, and then a reverse of the flow takes place as a result of the drum wall. The lowest particle velocity on the profile slightly depends on the drag model and[START_REF] Syamlal | The derivation of a drag coefficient formula from velocity-voidage correlations[END_REF]'s correlation seems to overestimate it in comparison with the other two drag models. Note that the used conditions ensure a rolling regime, which generates a tumbling motion of solids within the rotary drum. The active region occupies a large part ensuring a maximum exchange of intra-particles in the upper zone of the bed.The predictive capabilities of CFD models in rotating drums has been assessed in many experiments.Santos et al. (2013b) showed that the employed drag models do not impact the CFD model prediction and, therefore, can be neglected. Other studies[START_REF] Tamburini | Influence of drag and turbulence modelling on CFD predictions of solid liquid suspensions in stirred vessels[END_REF];[START_REF] Upadhyay | An assessment of drag models in Eulerian-Eulerian CFD simulation of gas-solid flow hydrodynamics in circulating fluidized bed riser[END_REF] highlight the crucial role of the drag model in the flow hydrodynamics

Figure 8 . 5 :

 85 Figure 8.5: Fertilizer particles velocity distributions predicted using the three drag models, for a rotational speed of 8 rpm and 30% filling degree.

  Fig. 8.6, the effect of the drag model on the evolution of the solid bed distribution is not significant. It is worth noting that an important phenomenon occurs at this stage: the upper percentage of fertilizer particles that reaches the top of the drum can not resist to ensure a flow continuity indicating an avalanche wave, and particles returns to the granular bed, which impacts directly the mixing quality within the drum and furthermore, favor inter-particles interactions for coalescence, breakage, and other mechanisms.

Figure 8 . 6 :

 86 Figure 8.6: Front view of a middle transversal plane of particles axial distribution with different drag models at different simulation time with 30% fill rate α s (particles volume fraction), a domain rotational speed of 8 rpm, and an inclination of 3.57 • .

Figure 8 . 7 :

 87 Figure 8.7: Particles velocity profile (d=2mm) for 8, 16, and 24 rpm rotational speeds with different drag models tested, with a representative snapshot of the depth of the granular bed δ (distance from the drum bottom to the interface between granules and air)

Figure 8 . 8 :

 88 Figure 8.8: Axial snapshots of granular temperature patterns in drums rotating at 8 rpm, 16 rpm, 24 rpm, and 32 rpm with the three studied drag models, at end time of 200 s.

Figure 8 . 9 :

 89 Figure8.9: Volume fraction of the the granular solid phase for fill level of 20%, 30%, and 40% for drum rotating at 8 rpm.

  Fig. 8.10 represents the fertilizer particles averaged axial velocity along the rotating drum for different operating conditions, with a constant fill rate of 30%. It is clear that the particles axial velocity increases gradually along the drum. As a result the axial component of the particles velocity is more and more important by mass conservation.On the other hand, the effect of the inclination degree on particles dynamics is varying along the drum length. By increasing the inclination slope, particles tend to occupy the lower end of the drum, due to gravitational force component in the axial direction. Fig.8.10 further displays that any increase of the inclination degree impacts the axial and radial particles distribution along the drum, which will enhance mass and heat transfer mechanisms between particles present in the active layer at the end part of the drum. The rest of particles defining the passive region can present a dead granular bed depending on the fill rate. Conjugated effects between inclination degree and rotational speed are clearly observed in Fig.8.10. By increasing the rotational speed and the drum inclination, the particles volume fraction decreases in the inlet zone.

Figure 8 . 10 :

 810 Figure 8.10: Axial particles velocity (m/s) distributions in divided zones per drum length (m) for different tested rotational speed(8 rpm, 16 rpm, and 24 rpm respectively)

  List of Figures 3.13 Relative power demand for a PBT inclined with 45 • from (Warmoeskerken and Smith, 1986) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 A schematic diagram Power number as a function of Reynolds number for the main turbine designs, from (Meyers, 2002) . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . 3.11 Turbulence Chemistry Interactions models matrix, from (Manual, 2009) . . . . . . 3.12 Flow map presentation for a Rushton turbine . . . . . . . . . . . . . . . . . . . . . xviii

  Axial snapshots of granular temperature patterns in drums rotating at 8 rpm, 16 rpm, 24 rpm, and 32 rpm with the three studied drag models, at end time of 200 s.8.9 Volume fraction of the the granular solid phase for fill level of 20%, 30%, and 40% for drum rotating at 8 rpm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	xx			List of Figures xxi xxv
	TVA S m ext		Tennessee Valley Authority Source term External	kg/m 3 .s
	UI p lf t		Uniformity Index Static pressure Lift	P a
	CM τ v m		Compartment Model Stress tensor Virtual mass	N /m 2
	EE J i gs		Eulerian-Eulerian Diffusion flux Gas-solid	Kg/(m 2 .s)
	EL D i,m f r		Eulerian-Lagrangian Diffusion coefficient Frictional	m 2 /s
	STCR µ bf RT STR µ t Sc t col kin CHAPTER 1	Stirred Tank Continuous Reactor Dynamic viscosity Unbaffled Rushton Turbine Stirred Tank Reactor Turbulent viscosity Schmidt Number Collisional Kinetic	kg / ms LIST OF TABLES N.s/m 2 -
	SIMPLE ρ unbf		Semi Implicit Method for Pressure Linked Equations Density kg / m 3 Baffled
	RTD k B		Residence Time Distribution Turbulent Kinetic Energy Bulk	m 2 /s 2
	RPM ε T		Rotation Per Minute Turbulence dissipation rate Turbine	m 2 /s 3
	PDF ω I		Probability Density Function Specific turbulent dissipation rate Impeller	m 2 /s 3
	tet G k		Tetrahedral production of turbulent kinetic energy	m 2 /s 2
	hex f k		Hexahedral Interaction forces between phases	N
	FVM p A		Finite Volume Method Partial pressure a component A	P a
	RST H A		Reynolds Stress Tensor Henry constant for a component A	-
	TPDF C A		Transported Probability Density Function Concentration of a component A	mol/l
	LDV C * A		Laser Dopper Velocity Equilibrium concentration in the liquid	mol/l
	k L		Mass transfer coefficient	m/s
	N		Impeller rotational speed	rpm
	n(d, z, t) Roman & Greec Symbols Number of density function P Power consumption B a Particle birth t Time D a Particle death u i Velocity component on i direction G Layering term Y α Masse fraction of spice α C d Drag coefficient ε g Gas holdup y + Structure of the boundary layer	--W -s -m/s ----
	K l		Mass transfer coefficient	h -1
	E Subscripts	Young's modulus	N /m 2
	F D		Force Drag	N
	g f		Gravity Fluid	m/s 2
	I f p		Moment of inertia Fluid-particle interaction	J/kg
	K i l		Reaction rate constant for specie i Liquid	1/s
	l g		Length Gas	m
	m p		Mass Particle	kg
	ṁ s		Mass flow rate Solid	kg/s
	8.10 Axial particles velocity (m/s) distributions in divided zones per drum length (m) Co Courant number -i Specie i
	N R	for different tested rotational speed(8 rpm, 16 rpm, and 24 rpm respectively) . . . Rotational impeller speed rpm Radius
	d t		Impeller diameter Tangential	m
	H 0		Reactor height Initial value	m
	N p b		Power number Bubble	-
	N q v		Pumping number Vapor	-
	T p		Impeller revolution period Power	s
	v m		Velocity Mixing	m/s
			xix xx xxiv xxv xxvi

• . 8.7 Particles velocity profile (d=2mm) for 8, 16, and 24 rpm rotational speeds with different drag models tested, with a representative snapshot of the depth of the granular bed δ (distance from the drum bottom to the interface between granules and air) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.8 5.1 The preneutraliser's operating parameters . . . . . . . . . . . . . . . . . . . . . . . 5.2 Solver setup for the case of preneutralizer simulation . . . . . . . . . . . . . . . . . 5.3 Details of torques calculated using CFD results . . . . . . . . . . . . . . . . . . . . 5.4 Details of consumed power using CFD results . . . . . . . . . . . . . . . . . . . . . 6.1 Parameters used for the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Details of the different grids used for grid independence study . . . . . . . . . . . 6.3 Hydrodynamical characteristics of the flow for 68 rpm, and 120 rpm . . . . . . . . 6.4 Parameters of the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Physical properties of the multiphase model . . . . . . . . . . . . . . . . . . . . . .

  Shear rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.6.6 Mixing time Θ t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.6.7 Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

	THEORETICAL ASPECTS OF STIRRED REACTORS:
	HYDRODYNAMICS, MASS TRANSFER AND NUMERICAL
	MODELING
	This chapter presents a literature review in the field of multiphase flows and mass transfer
	phenomena to illuminate the state of the art of preneutralization industrial unit; that is a gas-liquid
	sparged chemical stirred reactor. It also presents a survey on the existing technologies for the
	preneutralization, details of the numerical setup of CFD models and the available advanced CFD
	methods for reacting multiphase flows modeling and simulation.
	Contents
	3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
	3.2 Theoretical aspects on industrial mixing . . . . . . . . . . . . . . . . . . . . . . 24
	3.3 CFD approach for complex flows modeling . . . . . . . . . . . . . . . . . . . . . 25
	3.3.1 Impeller rotational movement modeling . . . . . . . . . . . . . . . . . . . 27
	3.3.2 Governing equations and mathematical models . . . . . . . . . . . . . . . 29
	3.3.3 Meshing of solution domain . . . . . . . . . . . . . . . . . . . . . . . . . . 30
	3.3.4 Numerical resolution procedure . . . . . . . . . . . . . . . . . . . . . . . . 31
	3.3.5 Turbulence modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
	3.3.6 Multiphase flow model formulation . . . . . . . . . . . . . . . . . . . . . 36
	3.4 Mass transfer in stirred sparged reacting flows . . . . . . . . . . . . . . . . . . 38
	3.4.1 Overview on mass transfer in gas-liquid flows . . . . . . . . . . . . . . . . 38
	3.4.2 Bubbles aspects and distribution in liquid sparged gas stirred reactors . . 40
	3.5 Turbulence-Chemistry Interactions (TCI) for multiphase reacting flows . . . 42
	3.5.1 The Eddy-Dissipation Model . . . . . . . . . . . . . . . . . . . . . . . . . . 43
	3.6 Hydrodynamics of stirred sparged gas liquid reactors . . . . . . . . . . . . . . 44
	3.6.1 Reynolds number Re & Froud number F r . . . . . . . . . . . . . . . . . . 45
	3.6.2 Power consumption P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
	3.6.3 Power number N p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
	3.6.4 Pumping capacity Q p & pumping flow number N q . . . . . . . . . . . . . 47
	3.6.5

  and r max being the larger of r i and r j . Variables a, b and c depends to the case study. many correlations exists as function of r max and V bond The viscous force consists of a normal and tangential component and is described by Nase et al.(Nase et al., 2001) as

					-δ n r max	+ b n	(4.16)
	with liquid surface tension σ F vis = 6πµR * R * -δ n	δn +	8 15	ln	R * -δ n	+ 0.9588 δt ,	(4.17)

.15) with capillary force F cap and viscous force F vis . The capillary force acts in normal direction and is defined by Soulie et al. as F cap = πσ r i r j c + exp a

Table 5 . 1 :

 51 The preneutraliser's operating parameters

	Parameter	Value	Unit
	Ammonium acid density	0.6894	Kg/m 3
	Phosphoric acid density	1272.64	Kg/m 3
	Ammonium acid viscosity 1.015 e -5 Kg/m.s
	Phosphoric acid viscosity	0.0203	Kg/m.s
	Rotating speed of agitator	68	rpm

table 5 .

 5 Table 5.2: Solver setup for the case of preneutralizer simulation

	Turbulence Model	k-ϵ model
	Fluids in the stirred vessel (Preneutralizer) Phosphoric Acid & Ammonia
	ézq Boundary conditions	Rotation MRF (68 rpm)
	Discritization-pressure	PRESTO!
	Discritization-momentum	Second order upwind
	Discritization	Species transport approach
	Discritization-turbulent kinetic energy	Second order upwind
	Discritization-turbulent dissipation rate	Second order upwind
	Residuals	10e -4

2:

ANSYS FLUENT can model the mixing and transport of chemical species by solving conservation equations describing convection, diffusion, and reaction sources for each component species. Multiple simultaneous chemical reactions can be modeled, with reactions occurring in the bulk phase (volumetric reactions) and/or on wall or particle surfaces. Species transport

Table 5

 5

		.3: Details of torques calculated using CFD results
	Designation	Torque (N.m) of unbaffled PN Torque (N.m) of the baffled PN
	Agitator	14501	27342.1102
	Shaft in static zone 0.011975	0.018656
	Shaft in MRF zone 0.000153	0.015900
	Total	14501.012128	27342.120246

Table 5 . 4 :

 54 Details of consumed power using CFD results

	Designation	Power (Kw) of the unbaffled PN Power (Kw) of the baffled PN
	Agitator	103.208	194.123
	Shaft in static zone 8.523006e -5	13.201e -5
	Shaft in MRF zone 1.088e -6	1.129e -3
	Total	103.208	194.466

Table 6

 6

	.1: Parameters used for the study
	Simulation parameters	Symbol Values	Unit
	Height of the reactor	H	7.3	m
	Width of the reactor	G	3.65	m
	NH 3 Nozzles diameter	D	0.105	m
	Inlet velocity at the central v in Ammonia gas density ρ	3.55 1.2	m/s kg/m 3
	Gas dynamic viscosity	µ	1.810 -5 Pa • s
	ACP density	ρ	2610	kg/m 3

Table 6 . 2 :

 62 Details of the different grids used for grid independence study Grid Resolution Number of Number of min / max of y + time (hr)

			nodes	cells	skewness
	1	Coarse	1066150	775058	0.44/ 0.37	49 120
	2	Medium	1599195	1162588	0.35/ 0.23	37 210
	3	Fine	1760629	1273690	0.30/ 0.16	32 270

Table 6

 6

	.3: Hydrodynamical characteristics of the flow for 68 rpm, and 120 rpm
			68rpm	120rpm
	Reynolds number (10 3 ) 3007,60	4420,70
	Power	(W) 21612,07 35042,41
	Power number	-	0,1	0,051
	Pumping number	-	0,69	0,55

Figure 6.6: Mass transfer profile in different levels (Z=3.7m; Z=2.7m; Z=1.4m) for 68 rpm and 120 rpm

Table 6 . 4 :

 64 Parameters of the study

	Simulation parameters	Symbol Value Unit
	Height of the reactor	H	7.3	m
	Width of the reactor	W	3.65	m
	NH 3 injection area diameter d NH 3 gas density ρ g NH 3 dynamic viscosity µ g ACP liquid density ρ l	0.104 m 1.215 kg/m 3 1.8e -5 P a.S 2650 kg/m 3
	ACP liquid viscosity Operating temperature	µ l T	0.021 P a.S 80 • C

  gs is the energy exchange between the fluid and solid phases, k θ s is the diffusion coefficient, and γ θs is the collisional dissipation energy given by[START_REF] Lun | Kinetic theory for granular flow of dense, slightly inelastic, slightly rough spheres[END_REF]:

	γ θ s =	12 1 -e 2 ss α 2 s ϱ s g 0 θ 3/2 s d s √ π	,	(8.14)
	and
		θ s =	1 3	u 2 s,i ,	(8.15)
	where e ss is the restitution coefficient for particle collisions, taken as 0.95, and the radial distri-
	bution function g 0 represents the probability of interparticle contact. Neglecting the convective
	and diffusive terms in Eq.(8.13) allows to reduce the transport equation for the granular tem-
	perature to the algebraic form
	(-p

.13) where (-p s I + τ s ) : ∇v s represents the resultant energy from the solid stress tensor, ∇ • (k θs ∇θ s ) is the diffusion of energy, Φ s I + τ s ) : ∇v sγ θ s -3K sg Φ gs = 0.

(8.16)

Table 8 . 1 :

 81 Physical properties of simulations

	Physical property	Value
	Air density (ϱ g ) Air viscosity (µ g )	1.255 kg/m 3 1.7894 × 10 -5 kg/m•s
	Fertilizer particles density (ϱ s )

Table 8 . 2 :

 82 Details of CFD simulated scenarios

	Parameter		Tested options
	Rotational speed (rpm)	8	16	24
	Inclination degree ( • )	3.5	5.5	7.5
	Filling degree (%)	20	30	40
	Drag model	Gidaspow (1994) Syamlal and O'Brien (1987)

dustries chimiques, en occurence dans la pétrochimie, pour la formation des opérateurs, la conception et l'optimisation des procédés. Cependant, il y a un manque des travaux de modélisation et de simulations numériques rigoureuses dans l'industrie des engrais phosphatés. Il existe de nombreux défis dans les systèmes de production d'engrais phosphatés, notamment les systèmes multiphasiques, telles que les liquides, les solides et le gaz, les particules de différentes tailles et distributions, et la variation des propriétés physiques, y compris la rhéologie et la thermodynamique. Malgré les ressources disponibles, peu d'études ont été consacrées à la compréhension de l'hydrodynamique au sein des lignes de production des engrais, afin d'augmenter le rendement de la production et d'adapter les conditions opératoires du procédé, de production à la qualité des matières premières. La dynamique des fluides numérique (CFD) représente un outil incountournable pour valider des modèles rigoureux et permettre une compréhension exhaustive de la sensibilité du procédé aux conditions de fonctionnement, et des effets des détails géométriques sur l'hydrodynamique de l'écoulement. De plus, la CFD permet une meilleure compréhension des phénomènes physiques, des mécanismes d'échange qui ont lieu dans les procédés étudiés, ce qui est casément pas possible par des expériences et dans les installations industrielles réelles ou pilote. C'est un outil efficace pour le scale up et pour l'optimisation des conditions de fonctionnement des unités chimiques. L'objectif de cette thèse est d'étudier la dynamique des fluides dans des pièces maîtresses du procédé de production des engrais Di-Ammonium Phosphate: la preneutralisation et la granulation, en utilisant le logiciel commercial de simulation Ansys Fluent, en identifiant l'effet des phénomènes complexes dans la formation d'hétérogénéités des flux multiphasiques, et en évaluant l'effet des conditions opératoires sur l'hydrodynamique des écoulements. De plus, à la lumière des résultats obtenus, optimiser la conception et le contrôle des unités de preneutralisation et de granulation et proposer de nouvelles recommandations pour une meilleure

General Introduction

Fertilizer chemical manufacturing process

Theoretical aspects of stirred reactors: hydrodynamics, mass transfer and numerical modeling

Theoretical aspects of rotating drums: hydrodynamics, and numerical modeling

Baffles effects on the hydrodynamics within the preneutralizer using CFD approach

CFD investigation of hydrodynamics in the preneutralizer : Analysis of the reacting flow reactor for DAP production

In silico CFD investigation of the granulation hydrodynamics in rotating drum

Conclusions

ACKNOWLEDGEMENTS

Firstly, I would like to thank Almighty ALLAH for his generous blessing upon me during all my educational journey.

My first thanks and deepst gratitude go to Prof. My Saddik Kadiri, Prof. Khalid Achoubir and Prof. M'barek Bakkas for their acceptance of reporting that work, and for taking part in my

Resources for the numerical calculations were provided by SIMLAB High Performance Calculations center, and Toubkal Super-calculator parts of Mohammed VI Polytechnic University HPC resources. I would like to thank IT Department and staff of MSDA program for their help with the procedures and set-up. The work leading to this PhD thesis has received funding from the OCP Group in partnership with Mohammed VI Polytechnic University

CHAPTER 4

THEORETICAL ASPECTS OF ROTATING DRUMS:

HYDRODYNAMICS, AND NUMERICAL MODELING

The aim of that chapter is to elaborate a survey on the theoretical aspects of granulation mechanisms, and the existing numerical models that describes the granulation and their resolution methods. On the other hand we presents the industrial equipments and technologies used for the granulation, With a special and detailed description of the main numerical methods for multiphase flow modeling with the presence of a granular phase. Also, the flow regimes classification theory is presented as the key parameter for the hydrodynamics analysis of dense flows using CFD modeling.

CHAPTER 5

BAFFLES EFFECTS ON THE HYDRODYNAMICS WITHIN THE PRENEUTRALIZER USING CFD APPROACH

This chapter presents a tailored 3D multiphase CFD model of the simplified preneutralizer chemical reactor. A benchmark study between the baffled and the unbaffled preneutralizer geometries will be presented to evaluate the hydrodynamics and its sensitivity to the CAD design. Especially, the vortex appearance dependence to the baffles consideration in the simulation domain. where ⃗ F q , ⃗ F lif t t, s, and ⃗ F vm are the external body force, lift force and virtual mass force, respectively. Only drag and gravity are considered in the present investigation. Details on the air and solid momentum exchange coefficient, the particles phase stress description, the collision and drag forces, also the particles characteristics integrated by the KTGF are given following the model of Santos et al. (2013b).

Numerical simulation conditions

The CFD code used for this study is the commercial ANSYS Fluent 2020 R2, and calculations have been developed to predict the flow behavior of a multiphase bed composed by granules and air within the flow domain defined as a rotating cylinder representing the rotary drum granulator. The rotating drum is inclined by 3.57°, and with real industrial dimensions as defined in Table 7.1. At the initial state, the granulator (Gr) contains a bed of granules with a height of 1250 mm and air. The granules are supposed spherical and the initial volume fraction is 1 over the zone of the granular bed. The PSD is defined in the solid region by patching two superposed zones, each zone contains a granules bed with the unified diameter (Table 7.1). We simulated the system with three different rotational speeds: 8 rpm (rotation per minute), 16

rpm and 24 rpm. The physical characteristics of the studied system and the simulated phases are presented in the following table : The finite volume discretization, and a pressure-velocity Diameter of second granules bed 1.00e -4

using SIMPLE Algorithm were adopted and the coefficient of restitution (e ss ) used was 0.9. We choose a time step of 1.00e -5 s, and each simulation was performed for a period of 100 s; because it was found that in this time period there was no holdup variation between two rotations of the drum. We set the interpolation cell-centre values to the cell boundaries using the First-Order-Upwind-method. The Multiple Reference Frame (MRF) approach was used to define the 7.3. Results and analysis

Flow behavior and regimes transition for air-granules simulation

The rotary granulator is set at different rotational speeds (8, 16 and 24 rpm) in order to analyze the different regimes of the granular bed motion with a homogeneous particles distribution.

Results in Figure 7.3 show the transition between the known typical transverse bed motions:

rolling, cascading, and cataracting (Mellmann ( 2001)) during rotation.

Transverse flow behavior with particles size distribution (PSD)

The fertilizer granules bed was set by defining two zones with a unified diameter for each one.

The diameters of particles are respectively: 1.00e -3 mm , and 1.00e -2 mm. Figure 7.4 shows the initial state of two zones with solid fraction KTGF approach. Regimes transitions is studied within the flow behavior in an industrial-size granulator, and have been investigated using a 3D CFD model. In this chapter, the simulation results allows a prediction of the multiphase flow behavior with identification of the inclination angle and rotational speed as key parameters that control the mixing and the size enlargement of particles. The simulation results agrree well with the experimental results already done using by [START_REF] Rong | Numerical analysis of sizeinduced particle segregation in rotating drums based on eulerian continuum approach[END_REF]). This work is a the preliminar step to identify the flow hydrodynamics sensitivity to the process parameters. In the next chapter, a further advanced study of the sensitivity to the fill level also, and the drag models in the CFD model construction well be investigated.

CHAPTER 8

CONCLUSIONS

In this dissertation, multiphase reacting flows in the preneutralizer stirred reactor and rotary drum reactor were investigated numerically. In this regard, some light has been shed on hydrodynamic parameters, which are of great importance in terms of the design and enhancing the yield of the DAP manufacturing process through an optimal monitoring of the studied masterpieces.

Summary of the thesis

In this thesis, we worked on the two principal unit operations in the industrial manufacturing process of Di-Ammonium Phosphate fertilizers; that are: The preneutralization and the granulation in a rotary drum granulator. The motivation of our study was to understand and master the physicochemical properties of the DAP fertilizer slurry by enhancing the yield of the preneutralization and mastering the granulation to upgrade the quality of the produced commercial fertilizer granules.

During this work, in the first hand; we treated the preneutralization, that consists on mixing the ACP liquid and the sparged Ammonia gas in a stirred chemical reactor, with a special geometry. We developed gradually a complete workbench for a three dimensional CFD model, subcomponent. We crossed both energy consumption and process effectiveness to determine the convenient parameters of operating for two processes in phosphate industry.

A new ring sparger is designed inspired from bioreactors, and its effect is investigated on the preneutralizer performance, especially for mass transfer and gas holdup. The model provides good agreement with previous literature studies, and confirms our design recommendations.

The new designed sparger proves its efficiency in enhancing the ammonia dispersion within the preneutralizer. The mass transfer coefficient and gas holdup sensitivity to the sparger diameter was examined, and the actual diameter increasing leads to enhancing these key parameters, up to the ratio d sp /D = 0.8. Therefore, in the gas-liquid modeling, Eulerian-Eulerian multiphase approach was used, in which gas and liquid phases are considered as interpenetrating continua.

The established works open numerous perspectives to be studied in these processes and also other phosphate processes.

Secondly, for the rotary drum granulator, we worked on batch three dimensional inclined rotating drum reactor, with regards to the real industrial dimensions. The granulation essentially aims at forming agglomerates of a well-defined size that is generally larger than the initial particle size. Depending on the operating parameters, the physical properties of materials, and the geometrical designs, rotating drums showed different flow regimes for the particles dynamics.

To investigate the fertilizer particles dynamics using the Eulerian approach, different drag Its numerical behavior is similar to that of the k-ϵ models.

DNS is conceptually the simplest approach to the problem of turbulence.

An initial velocity field is assumed and iterated and the velocity field is evolved over time. DNS requires considerable resources.The scale for the DNS resolution of small eddies is at least 1000 times smaller than those used in LES. DNS is the leading edge of turbulence research. So far only basic problems are being addressed. Applications of DNS can be used to study fundamental properties of turbulence, provide data to compare simplified turbulence models, and simulate other effects non-Newtonian phenomena.