CFD-DEM Modeling and Simulation Coupled to a Global Thermodynamic Analysis Methodology for Evaluating Energy Performance: Biofertilizer Industry

: This work develops a methodology based on real chemical plant data collected from a Nitrogen-Phosphorus-Potassium fertilizer (NPK) cooling rotary drum. By blending thermodynamic variables given by global energy and mass balances with computational ﬂuid dynamics-discrete element method (CFD-DEM) modeling and simulation, the methodology provides an initial approximation to the understanding of heat transfer inside industry rotary coolers. The NPK cooling processwas modeledin CFDsoftware SimcenterSTAR − CCM + 13.06.011 usinga Eulerian–Lagrangian scheme through a coupled CFD-DEM method using one-way coupling. The average temperature of the NPK particles was obtained as well as the average mass ﬂow of the particles dropping as the drum was rotating. The analysis was performed for two-particle diameters (8 and 20 mm) during 17.5 s. The average heat transfer coe ﬃ cient between the ﬂuid and the NPK particles during the simulated time was obtained. A thermodynamic analysis was carried out using instantaneous energy and mass balances. Prandtl, Nusselt, and Reynolds numbers were obtained for each simulated time step. Finally, through a non-linear regression using the Marquardt method, a correlation between Prandtl, Nusselt, and Reynolds number was developed that allowed analyzing the rotating drum. Results showed that the proposed methodology could serve as a useful tool during the design and analysis of any given rotary cooler, allowing calculation of the heat transfer coe ﬃ cient and obtaining the process variables that could expand the machine operational capabilities due to the knowledge of the Nusselt number as a function of the drum working parameters.


Introduction
Countercurrent rotary drums are commonly used in the bio-organic fertilizer industry for cooling of granular solids that have high internal moisture retention, low thermal conductivity, and hygroscopicity, such as the case of organic fertilizers. They are often made of a long cylindrical shell that rotates upon bearings along its longitudinal axis, which is inclined to the horizontal to induce product motion from the inlet towards the outlet of the cylinder (Figure 1a) [1]. A blower is placed above the product outlet to provide the cold countercurrent airflow needed for cooling. Besides, to promote air-material contact, most drums have lifters, which are fin-like structures placed along the cylinder length that lift the material from the bed and showers it through the air stream as the drum rotates, creating a cross-flow heat transfer scheme (Figure 1a) [1].
The design of a rotary drum relies on designer expertise and experience with the product and the process and the application of certain theoretical principles and empirical correction factors [2]. Drum flights and shells design have an essential effect on the cooling process [3]. During operation, some particles are sliding and rolling along with the drum, while others are being lifted or falling in spreading cascades through the air stream and re-entering the bed at the bottom, giving rise to complex particle dynamics which affect fluid-solid interactions (Figure 1b) [4]. Nevertheless, a better understanding of the heat transfer phenomena in these processes is still needed for process optimization.
promote air-material contact, most drums have lifters, which are fin-like structures placed along the cylinder length that lift the material from the bed and showers it through the air stream as the drum rotates, creating a cross-flow heat transfer scheme (Figure 1a) [1].
The design of a rotary drum relies on designer expertise and experience with the product and the process and the application of certain theoretical principles and empirical correction factors [2]. Drum flights and shells design have an essential effect on the cooling process [3]. During operation, some particles are sliding and rolling along with the drum, while others are being lifted or falling in spreading cascades through the air stream and re-entering the bed at the bottom, giving rise to complex particle dynamics which affect fluid-solid interactions (Figure 1b) [4]. Nevertheless, a better understanding of the heat transfer phenomena in these processes is still needed for process optimization. NPK grade 15-15-15 (Fertilizer with 15% Nitrogen, 15% Phosphorus, 15% Potassium) is a watersoluble granular fertilizer of 2.7 mm average particle diameter commonly produced in the Colombian biofertilizer industry. After drying, NPK enters a countercurrent rotary cooler at circa 90 °C and must be below 45 °C at the outlet to decrease the likelihood of compaction and agglomeration during curing and packaging, which negatively affects the quality of the product. Operators often achieve outlet temperatures by adjusting operating conditions such as reducing drum inclination angle, drum rotation, NPK feed rate, and fuel mass flow. Once an industrial rotary drum is manufactured and put into operation, however, further geometrical and operational modifications for process optimization are harder to incorporate, since they generally require long plant downtimes not previously accounted for in maintenance schedules, as well as additional design and manufacturing costs for geometry changes. Thus, these measures ensure a decrease in product throughput and an increase in process costs, which ultimately diminish overall company profits. For these reasons, the industry permanently demands alternate low-cost approaches for process optimization.
Computational techniques have been increasingly employed for understanding complex solidfluid interactions in multiphase flows. Among them, the discrete element method (DEM) has been widely accepted as an effective method in addressing complex inter particle phenomena in problems involving granular materials, such as granular flows and powder mechanics [5]. In recent years, coupled DEM and computational fluid dynamics (CFD) have been used for modeling the interaction between highly packed low-diameter granulated solids and gas/liquid fluids in fluidized beds, NPK grade 15-15-15 (Fertilizer with 15% Nitrogen, 15% Phosphorus, 15% Potassium) is a water-soluble granular fertilizer of 2.7 mm average particle diameter commonly produced in the Colombian biofertilizer industry. After drying, NPK enters a countercurrent rotary cooler at circa 90 • C and must be below 45 • C at the outlet to decrease the likelihood of compaction and agglomeration during curing and packaging, which negatively affects the quality of the product. Operators often achieve outlet temperatures by adjusting operating conditions such as reducing drum inclination angle, drum rotation, NPK feed rate, and fuel mass flow. Once an industrial rotary drum is manufactured and put into operation, however, further geometrical and operational modifications for process optimization are harder to incorporate, since they generally require long plant downtimes not previously accounted for in maintenance schedules, as well as additional design and manufacturing costs for geometry changes. Thus, these measures ensure a decrease in product throughput and an increase in process costs, which ultimately diminish overall company profits. For these reasons, the industry permanently demands alternate low-cost approaches for process optimization.
Computational techniques have been increasingly employed for understanding complex solid-fluid interactions in multiphase flows. Among them, the discrete element method (DEM) has been widely accepted as an effective method in addressing complex inter particle phenomena in problems involving granular materials, such as granular flows and powder mechanics [5]. In recent years, coupled DEM and computational fluid dynamics (CFD) have been used for modeling the interaction between highly packed low-diameter granulated solids and gas/liquid fluids in fluidized beds, tumbling mills, and particle impaction in water reservoirs, among others [6][7][8][9][10][11]. Nevertheless, real industrial processes such as NPK cooling in a rotary drum involve billions of particles which require vast computational power, which is unfeasible in most research centers [12]. Also, the effects of discontinuous or granulated materials on fluid properties (also known as two-way coupling) requires additional convective terms in the fluid's conservation of mass, momentum, and energy equations. Hence, simulations become much more computationally intensive [13].
To overcome the above difficulties, a methodology was developed that blends thermodynamic variables given by global energy and mass balances with CFD-DEM modeling and simulation that permits an initial approximation to the understanding of heat transfer inside industrial rotary coolers. Initially, the process is modeled as the cooling process of NPK 15-15-15 inside a 50 T h countercurrent rotary cooler located in a Colombian biofertilizer industry. The drum geometry and lifting flights were based on the industrial-scale design. The rotating drum was analyzed with a Eulerian-Lagrangian scheme using a CFD-DEM one-way coupled physics model corresponding to the airflow and the NPK granulated solid, respectively. Two particle sizes were analyzed: 20 mm and 8 mm diameter NPK spherical particles inside the rotating cooler in a 1 m and 0.1 m drum section length, respectively. The novelties of this work are presented below: • A CFD-DEM computational model for NPK in the rotary cooler was developed to study the influence of NPK particle diameter on its average temperature along the drum section length; 17.5 s of simulated cooling time was required to validate the results.

•
The process average heat transfer coefficient throughout the simulated time was obtained. • A thermodynamic model of the drum operation was developed using a heat exchanger model, which was fed by the previously found average heat transfer coefficient of the process.

•
Through a non-linear regression method, a correlation of dimensionless numbers was obtained which determines, in a given interval, the process operation characteristics.
CFD-DEM models have been used to study many different complex problems involving particle-fluid flow interaction [14,15] and have been found useful to model fluidized beds. The first study using this approach [16] simulated plug flow through horizontal pipes. Since then, many different problems have been solved with this method [17][18][19].
This paper offers a novel approach that mixes thermodynamic analysis (boundary conditions for the cooler), knowledge about the rotary system (geometry, rotational speed), as well as transport phenomena conditions (air velocity, particle size), to develop a methodology able to solve complex problems in the chemical industry. We expect this work to serve as a quick supporting tool for the design or modification of countercurrent rotary coolers in the biofertilizer industry, particularly in cases where new granulated materials are to be employed, or additional product requirements are placed which require changes in the process current operational conditions.

Materials and Methods
The NPK cooling process was simulated in CFD software Simcenter STAR − CCM + V.13.06.011 using a Eulerian-Lagrangian scheme through a coupled CFD-DEM method. Momentum, heat, and mass transfer, were exchanged in only one direction to reduce computational costs, also known as one-way coupling. In this method, only the continuous Eulerian phase (air) influences the solid particles. Hence, the effects of the granulated solid on the air, such as displacement, interphase momentum, mass, and heat transfer, are not initially considered in the methodology. The average temperature of NPK particles, as well as their average mass flow falling from the top of the rotating drum during 17.5 s for two-particle diameters (8 and 20 mm), is obtained. These results were then used to calculate the average heat transfer coefficient between the fluid and the NPK particles during the simulated time. Afterward, a thermodynamic analysis was performed using energy, and mass balances of the simulated section of the rotating cooler to obtain the heat transfer coefficient of the process. Also, the instantaneous Prandtl, Nusselt, and Reynolds numbers of the process for each simulated time step were calculated. Finally, through a non-linear regression using the Marquardt method, a correlation between Prandtl, Nusselt, and Reynolds numbers was obtained for the analyzed rotating cooler.

Particle Dynamics
Linear motion description for the NPK granular flow with spherical solid particles is modeled using the DEM since it extends the Lagrangian formulation to account for inter-particle interactions in the particle equations of motion, which are essential in highly loaded flows. The equation of conservation of linear momentum for a DEM particle of mass m p is given by Equation (1), where V p denotes the instantaneous particle velocity, F s is the resultant of the forces acting on the surface particle, and F b is the resultant of the body forces. These forces are also decomposed according to Equation (2), where F d is the drag force, F p is the pressure gradient force, F g is the gravity force, F c is the contact force from the DEM, and F MRF is the force produced by the rotating reference frame.
For the solid-fluid interactions, the resultant F s forces represent the momentum transfer from the continuous phase to the particle. The drag force is given by Equation (3), where ρ is the density of the continuous phase, A p is the projected area of the particle, V s is the particle slip velocity, and C d is the drag coefficient of the particle given by the Schiller-Naumann correlation, which is suitable for spherical solid particles. The pressure gradient force F p is defined according to Equation (4), where V p is the volume of the particle and ∇p static is the gradient of the static pressure in the continuous phase.
For the particle body forces, the gravity force is given by Equation (5), where g is the gravitational acceleration vector. The contact force F c represents inter-particle and particle-boundary interaction and is presented in Equation (6), where F cm is the contact force model. A modification of the linear spring contact model developed by Cundall and Strack was used [20].
The normal and tangential forces are defined by Equation (7), where K n is the normal spring constant, K t is the tangential spring constant, N n is the normal damping, ν n is the normal velocity component of the relative sphere surface velocity at the contact point, C f s is the static friction coefficient, and d n and d t are overlaps in the normal and tangential directions at contact points.
The normal and tangential spring stiffness is given by Equation (8). E eq , G eq , and R eq are the equivalent Young's modulus, shear modulus, and radius of the interacting particles, and they are Processes 2019, 7, 673 5 of 20 calculated according to Equation (9), where E A and E B are Young's modulus of the particles, v A , and v B their Poisson's ratios, and R A and R B their radii.
The normal and tangential damping is given by Equation (10), where N n damp and N t damp are the normal and tangential damping coefficient, and M eq is the equivalent particle mass. The normal and tangential coefficients and equivalent particle mass are given by Equation (11), where C nrest and C trest are the normal and tangential coefficients of restitution defined by the physical properties of the material, and M A and M B the mass of each colliding particle.
The force produced by the moving reference frame is given by Equation (12), where ω is the angular velocity vector of the rotating reference frame and r is the distance vector to the axis of rotation.
Rotational motion for DEM particles is described by orientations and, therefore, their angular momentum must also be conserved, and it is represented by Equation (13), where I p is the particle moment of inertia described by a second-order tensor, w p is the particle angular velocity, M b is the drag torque, that is, the moment that acts on the particle due to rotational drag, and M c is the total moment from contact forces.
The drag torque reduces the difference between a particle and the fluid in which it is immersed, and is defined by Equation (14), where C R is the rotational drag coefficient. Ω is the relative angular velocity of the particle to the fluid and is given by Equation (15), where v is the fluid velocity, and w p is the angular velocity of the particle. The rotational drag coefficient is defined by Equation (16) where Re R is the rotational Reynolds number defined by Equation (17).
Processes 2019, 7, 673 6 of 20 The total moment from contact forces is defined according to Equation (18), where r c is the position vector from the particle center of gravity to the contact point, and M cm is the moment that acts on the particle from rolling resistance. Rolling resistance was modeled with the proportional force method with a defined coefficient of rolling resistance µ r .
In order to study the heat and mass transfer process inside the cooling rotary drum, the solid-fluid heat transfer needs to be defined properly. The equation of conservation of mass of a material particle is given by Equation (19), where . m p is the rate of mass transfer to the particle. This term is zero since no evaporation occurs. The particle energy balance is shown in Equation (20), where Q rad represents heat transfer by radiation and Q s represents heat transfer by other sources, and both of them are negligible in this process.
Convective heat transfer Q t is calculated according to Equation (21), where h is the heat transfer coefficient, and A s is the particle surface area. The heat transfer coefficient was obtained from the particle Nusselt number presented in Equation (22), where k is the thermal conductivity of the fluid. Nusselt number was obtained using the Ranz-Marshall correlation defined by Equation (23), where Pr is the Prandtl number of the air.
Q t = hA s T − T p (21) The solid-solid heat transfer process is taken into consideration. When particles make contact with each other or with the wall, heat is transferred through conduction. In this work, the drum wall was considered adiabatic. Thus, conductive heat transfer was only considered between particles. Particle-particle heat transfer is given by Equation (24), where r c is the contact radius and k is the equivalent thermal conductivity of the two particles and T i , T j the temperatures of particle i and j. The equivalent thermal conductivity is calculated according to Equation (25), where k i and k j are the thermal conductivities for particles i and j.
The impact heat model was employed for calculating the heat production that results from friction and damping in DEM particles. This model has a linear formulation given by Equation (26), where c t and c n are the fractions of frictional and damping work that are converted to heat, f t and f n are the frictional and damping forces on the particle, and v t and v n are the relative tangential and normal impact velocities. q r = c t f t v t + c n f n v n (26) Air was modeled as a constant density gas using a Eulerian scheme. The equations of conservation of mass, momentum, and energy were solved with mesh motion, which provides an additional flux in the convective terms. This set of equations is defined by Equation (27), where ν g is the grid velocity in the reference frame and ν r is the relative velocity with respect to the reference frame, σ is the stress tensor, f b is the resultant body forces (gravity), E is the total energy per unit mass, q is the heat flux and S u and S e are mass and energy sources When the mesh is moving, cells shape and position change with time. Hence, an additional equation was solved to enforce space conservation, described by Equation (28).

Turbulence Model
Turbulence was modeled with the realizable k-model, which exhibits superior performance for flows involving rotation than the k-turbulence model. This model added two equations which solve the turbulent kinetic energy k e and turbulent energy dissipation rate [21], which are given by Equation (29): In these equations, S is the modulus of the mean rate-of-strain tensor, P k represents the generation of turbulence kinetic energy due to the mean velocity gradients, calculated in the same manner as standard k-model, and P b is the generation of turbulence kinetic energy due to buoyancy, calculated in the same way as standard k-model [21,22]. C 1 , η and S and µ are given Equation (30), while the turbulent viscosity is defined by Equation (31) and the terms required are defined in Equation (32). The final terms defining the turbulence model are presented in Equation (32), where Ω ij is the mean rate-of-rotation tensor viewed in a rotating reference frame with the angular velocity ω k , and the model constants A 0 and A s are given by Equation (33).

Heat-Transfer Model
The average NPK temperature is used for each time step to obtain the heat-transfer coefficient between air and solid particles. NPK average temperature differential is given by Equation (34), where dt is the simulation timestep, T t+dt is the average NPK temperature in time t + dt and T t is the average NPK temperature in time t.
Heat transfer between NPK and air is assumed to occur mostly during PK falling from the flights during drum rotation. Hence, NPK thermal energy is defined by Equation (35), where m m f is the average NPK mass that falls during one time step. It is assumed that this energy is given to the flowing air, increasing its temperature accordingly. Thus, the air temperature differential is defined by Equation (36), where m a is the air mass passing during one time step, and Cp a is the air-specific heat. The final expression for air temperature after one time step is given by Equation (37), where T ai is the air temperature at the inlet.
Qm t+dt = m m f Cp m dT t+dt (35) In order to obtain the average heat-transfer coefficient h, the process was modeled as a counter-current heat exchanger. This method is suitable in cases where flow rates and heat transfer area are constant. Thus, a generic heat exchanger was employed with two ends: A, at which the hot material enters, and hot air leaves the system, and B, where cold material leaves and cold air enters the system. The logarithmic mean temperature difference LMTD is defined as the difference between the hot and cold feeds at each end of the countercurrent heat exchanger and is given by Equation (38).
The LMTD was used to determine the heat transfer coefficient h for each timestep in the heat exchanger model, and it is calculated according to Equation (39), where dQm dt is the heat exchanged between NPK and the air during the defined timestep in the heat exchanger model, and A f p is the area Processes 2019, 7, 673 9 of 20 of heat exchange corresponding to the total area of the particles falling during one timestep. This term, the rate of heat transfer dQm dt during the time step dt is given by Equation (40).
The area of heat exchange A f p is given by Equation (41).
The average heat transfer coefficient for each simulation was found by averaging the heat transfer coefficients of each timestep during the entire simulated time.

Thermodynamic Model
A global energy balance was used to obtain a generalized thermodynamic model of the cooler. First, the heat was assumed to transfer from the NPK material to the air without loses. Hence, the energy balance is presented by Equation (42) The air outlet temperature was then obtained from Equation (44).
The heat-transfer coefficient from Section 2.2 was used for both particle diameters in a counter-current heat exchanger model of the entire NPK industrial cooler. The inlet-outlet temperatures were used to in the thermodynamic model to find the exchanged heat in the cooling process, which is calculated from Equation (45).
A he is the total area for heat exchange, and LTMD he is the logarithmic mean temperature difference of the heat exchanger. The total area of heat exchange is defined by Equation (46), where L is the rotating cooler length, and l is the length of the simulated drum model.
The logarithmic mean temperature difference is given by Equation (47).
Through trial and error, NPK outlet temperature that makes Q he = Q m was calculated. NPK and air inlet and outlet temperatures obtained from this model were validated with the ones given as boundary conditions in the real industrial rotary cooler.

Dimensionless Numbers
The Nusselt number and Reynolds number for the NPK material, as well as the Prandtl number for the air for each particle diameter during each time step, were obtained. Particle Nusselt is given by Equation (49).
The Prandtl number of the air was found by linear interpolation of published Prandtl numbers for given air temperatures. The average air temperature was defined according to Equation (50).
NPK Reynolds number was obtained from Equation (51), where v max is the maximum air velocity throughout the drum, and ν is the air kinematic viscosity at the average air temperature. Air kinematic viscosity was obtained by linear interpolation of published kinematic viscosities for given air temperatures. The maximum air velocity is given by Equation (52), where V ai is the average air inlet velocity and V f is the void factor, which is defined by Equation (53). V d is the drum volume, and V m is the volume occupied by the NPK. The drum volume is calculated according to Equation (54), where D is the drum diameter, and l is the drum length in the simulation.
The volume occupied by the NPK material is given by Equation (55), where N is the number of NPK particles in the model.

Non-Linear Regression for Dimensionless Numbers Correlation
The values of the dimensionless numbers obtained for each particle diameter at each time step were used to find a correlation of the form given by Equation (56), where a, b, and c are constants that relate the dimensionless numbers. The Marquardt non-linear regression method was used in the software Statgraphics Centurion V16.1 to obtain the values of these constants.

Rotary Cooler Characteristics
A fully operating Colombian industrial NPK rotary cooler was used as a case study for the proposed methodology. The cooler has a diameter of three meters, is 19 meters long, rotates at 5 rpm, and is tilted down 1.5 • to the horizontal. Four different flight designs are employed in the industrial drum. The flights are arranged in a defined order along the drum length to complete a total of 13 sections. Atmospheric air is used as cooling fluid and enters the drum using a blower which generates a mass flow of 25  Table 1.

Model Description
Considering that computational resources for simulating NPK cooling of such a large number of solid particles were unfeasible to acquire, the drum length was reduced, and particle diameter was increased to obtain a total mass that achieved less than a million particles but still was representative of the cooling process. Hence, two models were developed, one of a one-meter drum section length with 20 mm diameter NPK spherical particles, and one of a 0.1 m section length with 8 mm diameter NPK spherical particles, both tilted down 1.5 • to the horizontal.
For 20 mm particles, drum length was reduced considering fluid dynamic conditions, particularly Reynolds number and Pressure coefficient. Reynolds number remained constant and was independent from the length of the drum. Since Reynolds was turbulent, the hydrodynamic entry region was 10 times the drum diameter, resulting in 30 m. With this condition, fully developed turbulent flow could not be achieved. For pressure coefficient, the pressure drop is proportional to the length of the drum presenting a linear behavior. Since a developing region will dominate the phenomenon inside, the length of the drum to be analyzed was dropped to 1 m, which represented 5.2% of the total length, which is well inside the developing region of the process. For 8 mm particles, Fourier number was considered because of the transient situation to be analyzed during the cooling process. For both particle sizes, 20 and 8 mm, Fourier number was held at 4 × 10 −3 . This allowed the calculation of drum length for 8 mm particles' size.
For these lengths, masses of 1315.789 kg and 131.578 kg of NPK were obtained, formed by 299,164 and 467,443 spheres. Both models used one of the four given types of lifters, which were chosen since they allowed the highest particle dispersion angle (angle between the first and last particle showering from the top of the drum). The geometry of the 0.1 m long rotating cooler and its lifters is shown in Figure 2.
Air inlet and outlet conditions were defined as mass flow inlet and pressure outlet in the software. Values for mass flows, pressures, and temperatures were defined as equal to those currently employed in the industrial cooler. Nevertheless, solid NPK particles did not flow from defined inlet and outlet boundaries. Instead, a random particle injector, a tool inside STAR-CCM+, was used to place particles inside the drum's fluid region in a randomized fashion with a defined initial temperature equal to the NPK inlet temperature of the industrial cooler. Hence, even though no NPK inlet flow was used, NPK particles were initially set with the industrial cooler NPK inlet temperature.
Wall boundary conditions were chosen for the drum carcass and phase-impermeable conditions for the air inlet and outlet boundaries, so particles could not escape the drum during cooling once they have been randomly injected ( Figure 3). Thus, a constant number of particles were randomly positioned inside the drum at the start of the simulation and could not leave the drum. Particle dynamics were obtained by solving the equations found in Section 2.1.1 for the constant number of particles defined in Table 1.

Simulation
The simulation was run for 17.5 s for each model, which corresponded to 1.45 rotations of NPK particles inside the drum. Longer simulation times did not provide additional information regarding NPK temperature drop dynamics since no considerable differences in average heat transfer coefficients were found. Furthermore, the logarithmic mean temperature difference behavior is well defined at this point, allowing the calculations of the heat to be removed inside the drum during the cooling process. The average air inlet velocity was used for calculating the maximum air velocity . Average NPK falling mass flow was calculated using the track model in STAR-CCM+, which allowed the tracking of each moving particle throughout time. A plane cutting the drum in half in the y-axis was defined, and the tracked particles which cross it throughout a specified time frame

Simulation
The simulation was run for 17.5 s for each model, which corresponded to 1.45 rotations of NPK particles inside the drum. Longer simulation times did not provide additional information regarding NPK temperature drop dynamics since no considerable differences in average heat transfer coefficients were found. Furthermore, the logarithmic mean temperature difference behavior is well defined at this point, allowing the calculations of the heat to be removed inside the drum during the cooling process. The average air inlet velocity was used for calculating the maximum air velocity . Average NPK falling mass flow was calculated using the track model in STAR-CCM+, which allowed the tracking of each moving particle throughout time. A plane cutting the drum in half in the y-axis was defined, and the tracked particles which cross it throughout a specified time frame

Simulation
The simulation was run for 17.5 s for each model, which corresponded to 1.45 rotations of NPK particles inside the drum. Longer simulation times did not provide additional information regarding NPK temperature drop dynamics since no considerable differences in average heat transfer coefficients were found. Furthermore, the logarithmic mean temperature difference behavior is well defined at this point, allowing the calculations of the heat to be removed inside the drum during the cooling process. The average air inlet velocity V ai was used for calculating the maximum air velocity v max . Average NPK falling mass flow was calculated using the track model in STAR-CCM+, which allowed the tracking of each moving particle throughout time. A plane cutting the drum in half in the y-axis was defined, and the tracked particles which cross it throughout a specified time frame were plotted ( Figure 4). The data was then exported to a CSV file, where each row corresponded to each crossing particle and each column to the crossing time and z-position. Particle falling mass flow for each model is given by Equation (57), where ∆t is the time frame between the first and last particle that passes through the were plotted (Figure 4). The data was then exported to a CSV file, where each row corresponded to each crossing particle and each column to the crossing time and z-position. Particle falling mass flow for each model is given by Equation (57), where ∆ is the time frame between the first and last particle that passes through the defined plane section. Measurements for were performed in time frames positioned after one second of real simulated time since this time was enough for the initially randomly localized particles to position inside the flights of the rotating drum and produce a quasisteady falling mass flow.
= # ∆ (57) Figure 4. Particle tracking in the transverse plane section of model one. Each circle corresponds to one particle passing through the defined plane in a given time point.

Mesh Generation
The fluid domain was discretized using three models: surface remesher, trimmer, and prism layer mesher. The first one resulted in a good quality surface mesh by adequately maintaining the drum and flights geometry. Trimmer mesh was employed for generating a predominantly hexahedral mesh with minimal cell skewness and a low number of trimmed cells. This method was sufficient for producing a mesh with appropriate curvature and alignment with the z-direction corresponding to the drum length ( Figure 5). Prism layer mesher was used to model the inclusion of a prism layer along the borders of the fluid domain to capture boundary layer flow. Mesh properties are shown in Table 2.

Mesh Generation
The fluid domain was discretized using three models: surface remesher, trimmer, and prism layer mesher. The first one resulted in a good quality surface mesh by adequately maintaining the drum and flights geometry. Trimmer mesh was employed for generating a predominantly hexahedral mesh with minimal cell skewness and a low number of trimmed cells. This method was sufficient for producing a mesh with appropriate curvature and alignment with the z-direction corresponding to the drum length ( Figure 5). Prism layer mesher was used to model the inclusion of a prism layer along the borders of the fluid domain to capture boundary layer flow. Mesh properties are shown in Table 2.

Time and Mesh Independence Analysis
Three different mesh sizes were generated for each model. Each finer mesh had 1.5 times the number of cells of the coarser ones. Convergence criteria were set to 10 for momentum, energy,

Time and Mesh Independence Analysis
Three different mesh sizes were generated for each model. Each finer mesh had 1.5 times the number of cells of the coarser ones. Convergence criteria were set to 10 −4 for momentum, energy, turbulent kinetic energy and turbulent energy rate of dissipation. Initially, a time independence analysis was performed. For it, Mesh number 1 was chosen for each model ( Table 3). The time step was decreased to half for each simulation, beginning with 0.01 and 0.05 for models one and two, respectively. Average NPK temperature was recorded at time t = 3.5 s for each simulation. Timesteps of 0.005 and 0.01 were chosen for model one and two since they showed relative errors lower than 1%. A mesh independence analysis was then performed using the chosen timesteps (Table 4). Average NPK temperature was again captured at time t = 3.5 s for each simulation with each finer mesh, and the relative error was calculated. Mesh number 2 was chosen for each model since it had a relative error lower than 1%. Table 5 shows the chosen meshes and timesteps for each model.  Figure 6 shows NPK 8 mm particles cooling inside a 0.1 m drum section during 4.5 s. Particles are lifted by the drum flights which move due to the angular rotation of the drum and fall at a given mass flow rate which is constant in the central plane defined in Section 3.3 of this work. From the figure, it can be seen that particles located in the middle of the lower bulk of NPK material show the lowest temperature (around 71 • C in blue), while particles near the flights and drum surface exhibit the highest temperature. This circumstance can be explained by the fact that air-particle interactions taking place near the center of the rotating drum achieve higher heat transfer rates due to airflow being fastest in this area as in near the rotating drum edges. Besides, these particles seemed to experience little movement compared to those in other drum regions. Even though some of them appear to add filling mass to the lower moving flights, most of them are not capable of entering the flights and move as a bulk material to the left by the action of contact forces by particles already inside the flights. Therefore, once they start going upwards, they end up falling again to lower flights already being filled by falling particles from the top flights of the drum. This phenomenon is not desirable since it could cause high temperature drops in some particles, while others could not achieve the desired outlet temperature, hence altering the overall quality of the product at the outlet.      Table 6 shows the heat transfer variables calculated for obtaining the model's heat transfer coefficient. Air masses are different due to different time steps, which is also expected for falling material mass. However, if timestep and length are made equal to model one in model two, material falling mass remains almost identical in both: ( = 0.01; = 0.04(2)(10) = 0.8 = 0.82 = ). Hence, the falling behavior is essentially the same for both particle diameters.    Table 6 shows the heat transfer variables calculated for obtaining the model's heat transfer coefficient. Air masses are different due to different time steps, which is also expected for falling material mass. However, if timestep and length are made equal to model one in model two, material falling mass remains almost identical in both: ( = 0.01; = 0.04(2)(10) = 0.8 = 0.82 = ). Hence, the falling behavior is essentially the same for both particle diameters.  Table 6 shows the heat transfer variables calculated for obtaining the model's heat transfer coefficient. Air masses are different due to different time steps, which is also expected for falling material mass. However, if timestep dt and length l are made equal to model one in model two, material falling mass remains almost identical in both: (dt = 0.01; m aM1 = 0.04(2)(10) = 0.8 Kg =0.82 = m aM2 ). Hence, the falling behavior is essentially the same for both particle diameters. The heat transfer coefficient is significantly increased for lower diameter particles, due to a greater area of heat exchange when the particle diameter is made smaller.  Table 7 shows the comparison between the thermodynamic models obtained for each simulated model and that found for the industry cooler. The thermodynamic behavior of model one provides a close approximation to air and NPK temperature boundary conditions when model one is used in the calculations. The calculated area of heat exchange for model one shows that its cooling capabilities are higher than the one found on the industry cooler. Thus, heat losses should be taken into account, as the simulation is set with adiabatic wall boundary conditions.

Correlation
Equation (58) shows the correlation found using the Marquardt method. Table 8 shows the analysis of variance performed on the model. The p-value states that the model parameters can explain the variation seen in the data. Therefore, the model is suitable for reproducing the data obtained from the thermodynamic balance and the DEM-CFD simulation.   Figure 9 shows the Nusselt vs. Reynolds for a fixed value of Prandtl. Note that Nusselt and Reynolds grow as particle diameter is increased. (58) 99.8634 Figure 9 shows the Nusselt vs. Reynolds for a fixed value of Prandtl. Note that Nusselt and Reynolds grow as particle diameter is increased.

Conclusions
This work describes and employs a new methodology that blends CFD-DEM modeling and global energy and mass balances to obtain an initial approximation of the heat-transfer phenomena taking place inside an industry-scale 50 counter-current rotary drum during the cooling process of NPK 15-15-15 biofertilizer. The low computational effort needed for simulating the proposed CFD-DEM model and the ease of the methodology offers a fast, concise method to elucidate process heat transfer dynamics in the current operating point. Besides, by the calculation of process dimensionless numbers and correlations, heat-transfer coefficients and overall process operating characteristics can be predicted when process variables such as drum angular velocity, air velocity, material inlet-outlet temperatures, and the average particle diameter are changed inside the range of valid Reynolds numbers for the defined correlation. Hence, process optimization and design modifications can be achieved without long plant downtimes given by the execution of new experimental designs. In this way, overall company profits are not severely affected.
Consequently, results from our simulation could be useful in many different cases. One scenario could be the following: There is a current rotating cooler operating at set parameters and the company wants to change variables such as cooler rotation, tilt angle, feeding rate, and particle diameter. After

Conclusions
This work describes and employs a new methodology that blends CFD-DEM modeling and global energy and mass balances to obtain an initial approximation of the heat-transfer phenomena taking place inside an industry-scale 50 T h counter-current rotary drum during the cooling process of NPK 15-15-15 biofertilizer. The low computational effort needed for simulating the proposed CFD-DEM model and the ease of the methodology offers a fast, concise method to elucidate process heat transfer dynamics in the current operating point. Besides, by the calculation of process dimensionless numbers and correlations, heat-transfer coefficients and overall process operating characteristics can be predicted when process variables such as drum angular velocity, air velocity, material inlet-outlet temperatures, and the average particle diameter are changed inside the range of valid Reynolds numbers for the defined correlation. Hence, process optimization and design modifications can be achieved without long plant downtimes given by the execution of new experimental designs. In this way, overall company profits are not severely affected.
Consequently, results from our simulation could be useful in many different cases. One scenario could be the following: There is a current rotating cooler operating at set parameters and the company wants to change variables such as cooler rotation, tilt angle, feeding rate, and particle diameter. After applying the methodology, the user can predict how these variables should be reprogrammed for achieving any given set of outlet boundary conditions, such as outlet temperature and outlet mass flow of granulated solid after cooling.