...........................................................................................................................................

Meneses-Chacón, Jaramillo-Ibarra and Mas de les Valls / INGE CUC, vol. 15 no. 1, pp. 133-141, Enero - Junio, 2019

Numerical analysis of the thermal and fluid dynamic behaviour of the flue gases in a traditional furnace for panela production

Análisis numérico del comportamiento térmico y fluidodinámico de los gases de combustión en un horno tradicional para la producción de panela

DOI: https://doi.org/10.17981/ingecuc.15.1.2019.12

Research Article. Reception date: 17/09/2018. Acceptance date: 30/01/2019.

Edxon S. Meneses-Chacón

Universidad Industrial de Santander, Bucaramanga, (Colombia)

edxon.meneses@gmail.com

Julián E. Jaramillo-Ibarra

Universidad Industrial de Santander, Bucaramanga, (Colombia)

jejarami@uis.edu.co

Elisabet Mas de les Valls

Universitat Politècnica de Catalunya, Barcelona, (Spain)

elisabet.masdelesvalls@upc.edu

.

How to cite this paper:

Edxon S. Meneses-Chacón; Julián E. Jaramillo-Ibarra; Elisabet Mas de les Valls; “Numerical analysis of the thermal and fluid dynamic behaviour of the flue gases in a traditional furnace for panela production,” INGE CUC, vol. 15, no. 1, pp. 133-141, 2019. http://doi.org/10.17981/ingecuc.15.1.2019.12

.

Abstract

Introduction− Panela is a product derived from sugar cane that is prepared using a traditional burner designed especially for this purpose. According to studies found in the literature, it was identified that the thermal efficiency of panela burners is 30% on average.

Objective− The objective of this investigation is to contribute to the search for new alternatives for the improvement of the low efficiency present on these systems, mainly affecting the flue gases duct.

Methodology− The development of this study is as follows: first, a research of the radiation and optical thickness effect in a simplified furnace is carried out. Afterward, a series of simulations with modifications in the design of the flue gas duct for a real size furnace are analyzed.

Results− The results showed that the radiation effect must be considered and, even though the optical thickness is low, it has a relevant impact in the heat transfer process due to the high temperatures in the furnace. A chaotic movement of the gases implied more heat transferred to the heaters and high values of Nusselt with the addition of new elements in the duct were obtained.

Conclusions− Arrangement 1, provides the best results with a Nusselt and thermal efficiency increase. No significant differences between the DOM and the P-1 radiation were found.

Keywords− CFD; turbulent flow; radiation heat transfer; industrial furnace.

.

Resumen

Introducción− La panela es un producto derivado de la caña de azúcar. En su elaboración se utiliza una hornilla tradicional, diseñada especialmente para este propósito. Según estudios encontrados en la literatura, se ha identificado que la eficiencia térmica de las hornillas paneleras se estima en un 30% promedio.

Objetivo− Esta investigación tiene como objetivo contribuir en la búsqueda de nuevas soluciones para el mejoramiento del nivel de eficiencia, modificando principalmente el ducto de humos.

Metodología− El desarrollo de este estudio es el siguiente: primero, se realiza una investigación del efecto de la radiación y del espesor óptico en un horno simplificado. Posteriormente, se realiza una serie de simulaciones con modificaciones en el diseño del ducto de humos para un horno de tamaño real.

Resultados− Los resultados mostraron que se debe considerar el efecto de la radiación. Aunque el espesor óptico sea bajo, tiene un impacto relevante en el proceso de transferencia de calor debido a las altas temperaturas en el horno. Un movimiento caótico de los gases implicó más calor transferido a las pailas, y se obtuvieron altos valores de Nusselt con la adición de nuevos elementos en el conducto.

Conclusiones− El arreglo 1, proporciona los mejores resultados con un aumento de la eficiencia térmica y de Nusselt. No se encontraron diferencias significativas entre los modelos de radiación DOM y P-1.

Palabras clave− CFD; flujo turbulento; transferencia de calor por radiación; horno industrial.

I. Introduction

Panela is a product obtained from sugar cane, which is known under different names: Panela in Mexico, Colombia and some countries in South America, Jaggery in India, Desi in Pakistan, and Naam Taan Oi in Thailand [1]. This product has been an important sweetening agent because of its characteristic taste, nutritional value and medicinal properties [2]-[3]. Panela is produced in 25 countries with a world annual production of 11.05 Million Tons (MT). India is the major producer with 70 % of worldwide production, manufacturing 6 MT annually [1]. Colombia was the first Panela producer in America and the second worldwide with a production of 1.5 MT annually in 2015 [4].

Fig. 1. Three-pan Panela furnace.

Source: Authors.

To produce Panela an especial furnace has been traditionally used (see Fig. 1). This unit consists of a combustion area, a flue gas duct, boiling pans, and a chimney. The flue gases, which are mostly obtained from the bagasse combustion, flow through the interior of the duct transferring thermal energy to the pans. These pans are open heat exchangers whose main purpose is to evaporate the water from the sugar cane juice.

Fig. 2. Energy balance in the furnace.

Source: Authors.

The performance of these furnaces, in terms of thermal and financial sustainability, is being questioned. Certainly, the thermal efficiency (part of the energy supplied by the fuel that reaches the pans for the water evaporation, in Fig. 2) could get up to 30 % in average and as low as 14 % [5].

The energy supplied by the fuel is not totally used to carry out the evaporation process, some of this energy is lost by [6]:

• Incomplete combustion.

• The exhaust gases.

• Through the walls and unburnt fuel.

The incomplete combustion and the energy lost can lead to some complications, such as environmental damage, bigger fuel consumption, longer time of batches, increase of the production costs, and in some cases, poor product quality [7].

An increment in the thermal efficiency, as well as fuel self-sufficiency, could be achieved by the customization of the furnace. For this, it will be crucial to implement an accurate heat transfer study of the actual process. Recently, some technological advances have been done mostly in the combustion chamber, the chimney and the boiling pans; all these modifications, experimentally obtained, have been applied in systems already in operation [8], [9], [10].

Numerical techniques have been used for the estimation of the fluid flow behavior and the heat transfer process in a few furnaces.

In the work by Osorio et al. [11] was reported a CFD study of the thermal evaluation of a sugar cane-burner. The results showed that the computational model is quite stable. They were able to predict the heat transfer process without noticeable differences compared with experimental data. For the solution of turbulence, they employed the Re-Normalization Group (RNG) k – ε model [12] with a mesh of 43509 control volumes.

La Madrid [13] did a numerical simulation of heat transfer phenomena on the open heat exchangers, and the results were compared to field measured data. They applied the Wilcox revisited turbulence model [14] as well as the P1 model [15] for the radiation part.

The specialized literature reflects a low level of research focused on the flue gas duct. As an additional contribution to the above mentioned numerical studies, the present work highlights the importance of the numerical tools in the process of improving thermally agricultural systems, the importance of analyzing the radiation effect of the flue gases [15], the comparison of the Discrete Ordinates Method (DOM) [16] and P-1 [17] radiation model, the usage of an open source code, in this case OpenFOAM [18], and the evaluation of new configurations based on the traditional flue gas duct geometry.

This work is composed firstly by a brief description of the mathematical model. Then, a comparison and further study of the radiation numerical models in a simplified industrial domain was performed. After that, simulations of an actual size Panela furnace two modifications of the original duct are showed. Finally, some conclusions are given

II. Mathematical Model

For most engineering applications, such as the one analyzed in this work, it is not necessary to resolve the details of the turbulent fluctuations. Therefore, time-averaged Navier-Stokes technique (RANS), introduced by Reynolds [19] that allows the computation of the mean flow without calculating the full time-dependent flow, was used. In these turbulence models, the equations are time-averaged before they are solved. Hence, and assuming an incompressible flow, they can be written as:

(1)

(2)

(3)

In these sets of equations i, , , μ and r represent the averaged velocity, temperature, and pressure, the dynamic viscosity and the radiative heat flux, respectively.

A. Turbulence model

As equations (2) and (3) show, a model for the Reynolds (turbulent) stress tensor 'j 'j and turbulent heat flux 'j ' is needed (closure problem). These additional terms appear as a consequence of the non-linearity characteristics of the convective term [20].

1) Boussines hypothesis

The Boussines approach is employed to calculate the turbulent stresses in equation (2), by means of a Linear Eddy-Viscosity Model (LEVM). This expresses that the Reynolds stress tensor is linearly related to the mean strains as follows [19]:

(4)

where k is the turbulent kinetic energy, μt the turbulent eddy viscosity, δij the Kronecker delta and ρ the density of the fluid.

2) Simple gradient diffusion

In this, the turbulent scalar flux is contemplated to be aligned with the mean scalar gradient. The turbulent diffusion coefficient is generally considered proportional to the eddy viscosity according to the turbulent Prandtl number, which is a specific constant of the turbulence model [21].

(5)

The k – ω version of Wilcox has been demonstrated to be particularly accurate [22] and has therefore been selected for use in the present study. The equations that constitute the k – ω model are [20]:

(6)

(7)

With closure coefficients and auxiliary relations as follows:

Table 1. Constant Values For The Standard Model.

.	α	β	β *	σ	σ *	ε
	5/9	3/40	9/100	1/2	1/2	β * ωk

Source: Wilcox [23].

The eddy viscosity can be expressed for the k – ω model as [20]:

(8)

B. Radiation models

Since gas radiation is considered in these simulations, besides the convective and conductive terms, the radiative effect is also present in the energy equation (see Eq. 3) as a source term ∇ ∙ r . The divergence of radiative heat flux is given by:

(9)

Where r– and s– are the position and direction of radiation intensity I, κ the absorption coefficient, and Ib(r–) = κ(T(r)4)/π. To obtain the radiation intensity field and then the radiative heat flux, the general equation of transfer for an absorbing, scattering and emitting gray medium must be solved. This equation is known as the radiative transfer equation (RTE) and is represented as [15]:

(10)

in which σs is the scattering coefficient, Φ(s– ∙ s–') is the scattering phase function for the radiation from incoming direction s–' and confined within the solid angle dΩ' to scattered direction s– confined within the solid angle dΩ. In this study, isotropic scattering medium is considered, in which the phase function is equal to unity. Due to the dependence on three spatial coordinates, two local directions coordinate and wavelength, an analytical solution is nearly impractical for most engineering applications. Thus, equation (10) has to be solved numerically using radiation transport models for spatial and directional dependencies [24].

Two numerical models have been tested for the solution of the RTE. One of them is the DOM. This method replaces the RTE by a discrete set of N differential equations that describe the radiation intensity field along N discrete ordinates or directions. The other method is the P-1 model. It is based on the reduction of the integral terms of the RTE to differential terms via a finite set of moment equations. The P-1 radiation model is the simplest case of the P-N model. To develop the general P-N method, the intensity at each special position is expressed as an expansion in a series of orthogonal harmonics and the series after a finite number of N terms is truncated. The mathematical background of the RTE resolution for each case can be found in [15].

III. Influence of radiation model on heat transfer for simplified configurations

Gases emit and absorb radiation at several narrow wavelength bands. The emission and absorption characteristics of the composing gas mixture also depend on the pressure and temperature [17]. In order to keep complexities to a manageable level, in this study the gas is assumed to be grey, and a total absorptivity and emissivity determined by an averaging process are used. Charts for the total emissivities given by Hottel were employed because of its wide application in radiation calculations with good accuracy [25].

According to [13], the mass fractions of the flue gases product of the bagasse combustion are:

YCO = 0.0001

YCO 2 = 0.0509

YH 2O = 0.037

YO 2 = 0.1472

YN 2 = 0.7648

Some of these gases may not cooperate with the radiation energy exchange. Molecules such as N2 and O2 are substantially transparent to radiation. On the other hand, gases with molecules such as H2O and CO2, may participate in the radiation process by absorption at moderate temperatures, and by absorption and emission at high temperatures [25]. For that reason, the flue gases considered in this investigation are the H2O and the CO2. The effect of the CO is neglected due to its low concentration in the mixture. The optical thickness was estimated following the procedure proposed by [25]:

Table 2. Heat Transfer Radiation Conditions.

.	Parameter	Value
	L	3.6V/As = 0.46
	Operating pressure	P = 1 [atm]
	Absorption coefficient	kλ = 0.33
	Optical thickness	τ = kλ * L = 0.15
	Albedo coefficient	ωa = 0

Source: Authors.

As reported by [17], special care must be taken due to ability of the medium to participate in the exchange of thermal energy. Therefore, in order to analyse the radiation effect a simulation in an industrial-like domain with only one pan was completed. Additionally, a comparison of the DOM and P-1 radiation model is performed since the former is more expensive computationally speaking, and the latter is generally not accurate for optically thin mediums (τ < 1) according to [16]. The discrete ordinates method (DOM). The dimensions are shown in figure 3.

Three simulations were carried out:

• The first one is done without considering the radiation effect (results denoted as No Radiation in figure 4). The boundary conditions applied are those explained in section IV-A, tables 3 and 4.

• The second uses the same boundary conditions, but the effect of radiation is considered using the DOM approach and an optical thickness of 0.15. The emissivity values used are: pans (stainless steel) = 0.4 and walls (refractory brick) = 0.94. (referred as Radiation τ = 0.15 DOM, in Fig. 4).

• The last simulation makes use of the same conditions as the second one, but the P-1 radiation model was applied. (referred as Radiation τ = 0.15 P-1 in figure 3).

Fig. 3. Furnace with one heater (dimensions in mm).

Source: Authors.

At x = 3.55 m, temperature profiles were taken (Fig. 4). Significant differences at the exit of the duct are observed. The temperature decreases due to the optical thickness. The elevated temperature of the flue gases allows the radiation intensity to gain energy as the beam crosses through the gas. This earned energy is lost by the medium. It is evident as well that the difference between the results of both radiation models is negligible. For this reason, the P-1 model is used in the industrial configuration, since it requires significantly less computational effort (around 70 %) than the DOM approach.

Fig. 4. Temperature profiles: Duct exit, x = 3.55 m, y = 0 m.

Source: Authors.

IV. Numerical simulation of a 3D industrial furnace.

In this section, the authors aim to present the numerical results of the thermal and fluid dynamic flow of the gases in an industrial furnace employed for producing Panela. A sketch of the geometry selected is shown in Fig. 5.

Fig. 5. Isometric view of the Furnace.

Source: Authors.

A. Boundary conditions

In the simulation, inlet boundary conditions were selected in order to emulate a real furnace able to produce 50 kilograms of Panela per hour (Kg/h). The dimension and the boundary condition are based on the work by Gordillo and Garcia [26].

Table 3 and 4 show the boundary conditions (hereafter referred to as patches due to the keyword used in OpenFOAM) for the configuration selected (Fig. 5). No slip boundary condition in the solid patches’ heaters, upper walls, ground, and walls are imposed. For temperature, a linear variation of 771.7 °C to 458.8 °C along duct is established due to the heat losses that are presented in the furnace bottom and walls boundaries [13]. A characteristic length l = 0.07 * L (L the hydraulic diameter) and a turbulence intensity (I) of 5 % is assumed. In the upper walls, a condition accounting the conductive effect of the bricks and the convective behaviors of atmospheric air is stablished (U).

Table 3. Boundary Conditions: Velocity, Temperature, And Pressure.

.	Patch	Velocity [m/s]	Temperature [°C]	Pressure (P/p) [m2/s2]
	inlet	ux: 0, uy: 0, uz: 6	T: 970	∂P⁄∂xi = 0
	output	∂u–i ⁄∂xi = 0	∂T⁄∂xi = 0	P = 0
	heater	u–i = 0	T: 100	∂P⁄∂xi = 0
	upperWalls	u–i = 0	Global heat transfer Coefficient (U)	∂P⁄∂xi = 0
	ground	u–i = 0	Linear profile	∂P⁄∂xi = 0
	walls	u–i = 0	Linear profile	∂P⁄∂xi = 0

Source: Authors.

Table 4. Boundary Conditions: Turbulent Kinetic Energy And Specific Turbulence Dissipation. The Turbulence Intensity And L The Characteristic Length.

.	Patch	Turbulent kinetic energy [m2/s2]	Specific turbulence dissipation [1/s]
	inlet	k = 3/2(u–i I 2)	ω = √k/l
	output	∂k⁄∂xi = 0	∂ω⁄∂xi = 0
	heaters	Wall Function	Wall Function
	upperWalls	Wall Function	Wall Function
	ground	Wall Function	Wall Function
	walls	Wall Function	Wall Function

Source: Authors.

The boundary conditions for radiation are the same applied in the simplified furnace. Emissivity values of 0.4 (stainless steel) and 0.94 (refractory brick) with an optical thickness of 0.15 were established.

B. Mesh

Mesh generation is one of the most important steps to simulate CFD problems. Fundamentally, CFD simulations are affected by the quality of the mesh. For this case, an unstructured mesh is created using SALOME (version 7.7.1), with Netgen-3D settings Max. Size: 0.025, Min. Size: 0.0001, Growth Rate: 0.1. After performing a partition of the complete domain, the mesh is generated. Boundaries are named in SALOME in accordance with those in Tables 3 and 4. The mesh is exported from SALOME into UNV format. Finally, to create an OpenFOAM mesh the ideasUnvToFoam tool provided in OpenFOAM is used.

To get more accuracy and reduce the numerical errors, a grid independence study, as well as a refinement process was applied. A starting grid of 200.000 control volumes has been used and, by a grid-refinement process, the number of control volumes is augmented, the last mesh has been chosen of 1.400.000.

C. Numerical methods

The set of equations (explained in section two) are solved over the entire domain using a segregated PISO-like algorithm implemented in OpenFOAM [18]. Central differences are employed for the evaluation of diffusion terms. Whereas, convective terms are discretized by means of a second order Upwind scheme [27].

V. Results

Fig. 6, shows the streamlines of the flue gases in the duct coloured by velocity magnitude. It is noticed that immediately after the entrance of the furnace, a low-velocity recirculation at the bottom wall is generated. It is also observed that the gases flow around the external surfaces of the pans, with an increase of the velocity below the third one due to the area reduction. At this point, the maximum velocity is reached. Finally, the gases slow down because of the area expansion and a small recirculation zone downstream the third pan is created near the upper wall.

Fig. 6. Velocity streamlines [m/s].

Source: Authors.

Fig. 7, shows the distribution of the Nusselt number along the upper wall and the pans midplane.

It is observed that the Nusselt is slightly bigger in the third pan than in the other two. This is due to the high velocity and the turbulence degree of the gases around heater 3 (high values of turbulent kinetic energy).

Furthermore, it can be observed that there are some peaks of Nusselt number in every pan. In the first and third one, the highest values of Nusselt are located at the beginning while the second pan has it at the end. This can be correlated to the amount of thermal energy transferred in some specific areas of the pans, as can be observed in Fig. 8. The high concentration of heat is caused by the impinging of the gases on each heater as they flow through the duct.

Fig. 7. Nusselt number distribution along upper wall and pans. Plane y = 0.

Source: Authors.

Fig. 8. Heat flux on each pan [W/m2].

Source: Authors.

A. Analysis of alternative duct geometry.

Additional simulations were conducted modifying the duct geometry in order to enhance the energy efficiency obtained in the previous section.

Two different duct geometries were studied. The duct geometry is changed to achieve larger velocity gradients near the pans and with a greater turbulence level. The first one (figure 9a), hereafter referred as arrangement 1, consists of placing transversal plates at the bottom of the duct, perpendicular to the main flow direction of the flue gases. One of the plates is located between the 1st and 2nd pan, and the other one between 2nd and 3rd pan. As this design could face an accumulation of ashes and unburnt material between the plates, a new arrangement is studied (figure 9b). It comprises the insertion of a series of steps throughout the bottom wall (floor). The goal is the same as the first option, but it allows an easy way to clean the area.

Fig. 10, shows the streamlines colored with the velocity magnitude. More recirculation zones and stagnation of the flow in the zone below the 1st and 2nd heater is observed in both configurations.

(a)

(b)

Fig. 9. New arrangements: geometrical configurations.

Source: Authors.

(a)

(b)

Fig. 10. New arrangements: velocity streamlines [m/s].

Source: Authors.

Fig. 11, depicts the turbulent viscosity field. As desired, the flue gases flow nearer to the heaters with greater turbulent intensity. Higher values of turbulent energy in the pans 2 and 3 are observed for arrangement 1. The amount of turbulent viscosity in this arrangement is up to four times the viscosity of the base case in the same pans.

(a)
(b)	(c)

Fig. 11. Non-dimensional turbulent viscosity νT/ν in the midplane y = 0: (a) Base case, (b) Arrangement 1, (c) Arrangement 2.

Source: Authors.

From Fig. 12, it can be noticed that the temperature at the end of the duct is lower than in the base case, which primarily means that more energy is being transferred to the pans (although some heat is being lost through the upper walls, see figure 13 at x = 4 m).

(a)

(b)

Fig. 12. New arrangements: Temperature distribution [K].

Source: Authors.

The chaotic movement of the flue gases allows more contact and separation of the boundary layers that flow around pans 2 and 3. This increases heat transfer, which can be quantified in terms of the Nu. It can be seen in figure 14 that the Nu values in arrangement 1, is three times bigger in the second pan (interval x = 2.6 to x = 3.6) and twice the amount in the third one (interval x = 4.1 to x = 4.6).

Fig. 13. Nusselt number throughout the length of the duct (upper wall and pans) y = 0.

Source: Authors.

The energy efficiency obtained from this configuration is 13%. For arrangement 2, all pans receive more energy; however, in the second one, the difference is not as high as the first arrangement. Moreover, the transport of energy, in the first pan, was not increased in any of the arrangements. Based on the results of arrangement 1, a plate could be placed immediately after the entrance of the furnace, so that more turbulent flow is created, and more heat could be transmitted to this pan.

V. Conclusions

This paper presents a turbulent, 3D numerical investigation of a forced convection heat transfer and fluid flow in an industrial furnace for panela production. Different aspects related to new developments of the geometry of the duct were presented.

New arrangements were tested based on the results of the industrial case. Any form of gas acceleration was sought; the chaotic movement of the gases implied more heat transferred to the heaters. High values of Nusselt with the addition of new elements in the duct were obtained:

• Arrangement 1 provides the best results with a Nusselt and thermal efficiency increase. Additional improvements must be proposed if the stagnation of un-burnt material is presented.

• Arrangement 2 is a reasonable option. However, additional simulations employing few and higher steps are suggested. The implementation of this design could be more expensive than the first one for furnaces already in operation.

Radiation plays an important role in the heat transfer process due to the high temperature of the flue gases even though the optical thickness is very low. Significant changes were found considering the medium as participant. Comparison of the DOM and the P-1 radiation model for a simplified furnace domain was presented. No significant differences were found, although the P-1 model saves up to 70 % in computing time.

Funding

This article is derived from the project entitled “Evaluación térmica por métodos computacionales de un modelo de hornilla panelera”, funded by Universidad Industrial de Santander (UIS), Colombia. February 2015 – July 2018

References

[1] P. V. K. Jagannadha Rao, M. Das, and S. K. Das, “Changes in physical and thermo-physical properties of sugarcane, palmyra-palm and date-palm juices at different concentration of sugar,” J. Food Eng., vol. 90, no. 4, pp. 559–566, Feb. 2009. doi: https://doi.org/10.1016/j.jfoodeng.2008.07.024

[2] N. Singh, D. Kumar, S. Raisuddin, and A. P. Sahu, “Genotoxic effects of arsenic: prevention by functional food-jaggery.,” Cancer Lett., vol. 268, no. 2, pp. 325–30, Sep. 2008. doi: https://doi.org/10.1016/j.canlet.2008.04.011

[3] A. P. Sahu and B. N. Paul, “The role of dietary whole sugar-jaggery in prevention of respiratory toxicity of air toxics and in lung cancer,” Toxicol. Lett., vol. 95, Supplement 1, p. 154, Jul. 1998. doi: https://doi.org/10.1016/S0378-4274(98)80615-2

[4] H. García, A. Toscana, N. Santana, and O. Insuasty, Guía tecnológica para el manejo integral del sistema productivo de caña panelera. Bogotá, Colombia: Ministerio de Agricultura y Desarrollo Rural, Corpoica, 2007.

[5] K. S. S. Rao, A. Sampathrajan, and S. A. Ramjani, “Efficiency of traditional jaggery making furnace,” Madras Agric. J., vol. 90, no. 3, pp. 184–185, Jan. 2003. Available: http://www.panelamonitor.org/media/docrepo/document/files/efficiency-of-traditional-jaggery-making-furnace.pdf

[6] V. R. Sardeshpande, D. J. Shendage, and I. R. Pillai, “Thermal performance evaluation of a four pan jaggery processing furnace for improvement in energy utilization,” Energy, vol. 35, no. 12, pp. 4740–4747, Dec. 2010. doi: https://doi.org/10.1016/j.energy.2010.09.018

[7] K. González, Determinación de pérdidas energéticas y sus puntos críticos, en hornillas paneleras Ward-Cimpa en la hoya del río Suárez, Univ. Industrial de Santander, Colombia, 2010.

[8] P. Arya, U. K. Jaiswal, and S. Kumar, “Design based improvement in a three pan Jaggery making plant for rural India,” Int. J. Eng. Res., vol. 2, no.3, pp. 264-268, Jul. 2013.

[9] O. Mendieta, “Desarrollo de un modelo experimental para el coeficiente de transferencia de calor en el proceso de evaporación del jugo de caña de azúcar en un arreglo de película delgada,” Univ. Industrial de Santander, Bucaramanga, Colombia, 2012.

[10] G. B. Agalave, “Performance improvement of a single pan traditional Jaggery making furnace by using fins and baffle,” Int. J. Adv. Res. Sci. Eng., vol. 4, no. 4, pp. 85–89, Apr. 2015. Available: https://www.ijarse.com/images/fullpdf/1429353638_12_Research_Paper.pdf

[11] J. Osorio, H. Ciro, and A. Espinosa, “Evaluación Térmica y Validación de un Modelo por Métodos Computacionales para la Hornilla Panelera GP150,” Dyna, vol. 77, no. 162, pp. 237–247, Jun. 2010. Available: http://bdigital.unal.edu.co/5373/1/jairoosorio.2010.pdf

[12] D. Choudhury, Introduction to the renormalization group method and turbulence modeling, Lebanon NH, USA: Fluent Inc., 1973.

[13] R. La Madrid, D. Marcelo, E. M. Orbegoso, and R. Saavedra, “Heat transfer study on open heat exchangers used in jaggery production modules – Computational Fluid Dynamics simulation and field data assessment,” Energy Convers. Manag., vol. 125, pp. 107–120, Oct. 2016. Doi: https://doi.org/10.1016/j.enconman.2016.03.005

[14] D. Wilcox, “Formulation of the k-omega Turbulence Model Revisited,” in 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, Jan. 8–11, 2007. doi: https://doi.org/10.2514/6.2007-1408

[15] M. F. Modest, “Chapter 23 - Inverse Radiative Heat Transfer,” in Radiative Heat Transfer, 3rd Ed., pp. 779-802, Cambrigde, MA, USA: Academic Press, 2013. doi: https://doi.org/10.1016/B978-0-12-386944-9.50023-6

[16] G. Colomer, M. Costa, R. Cònsul, and A. Oliva, “Three-dimensional numerical simulation of convection and radiation in a differentially heated cavity using the discrete ordinates method,” Int. J. Heat Mass Transf., vol. 47, no. 2, pp. 257–269, Jan. 2004. doi: https://doi.org/10.1016/S0017-9310(03)00387-9

[17] M. F. Modest, Radiative Heat Transfer, 3rd Ed., Cambrigde, MA, USA: Academic Press, 2013. Doi: https://doi.org/10.1016/B978-0-12-386944-9.50023-6

[18] OpenCFD Ltd (ESI Group), “OpenFOAM.” .

[19] S. B. Pope, Turbulent flows. Cambridge, MA, USA: Cambridge Univ. Press, 2000. doi: https://doi.org/10.1017/CBO9780511840531

[20] D. Wilcox, Turbulence modeling for CFD. La Cañada, CA, USA: DCW Industries, Inc., 1998.

[21] J. E. Jaramillo Ibarra, “Suitability of different RANS models in the description of turbulent forced convection flows: application to air curtains,” TDX (Tesis Dr. en Xarxa), Univ. Politècnica de Catalunya. Dept. de Màquines i Motors Tèrmics Barcelona, España, 2008.

[22] J. E. Jaramillo, C. D. Pérez-Segarra, A. Oliva, and K. Claramunt, “Analysis of different RANS models applied to turbulent forced convection,” Int. J. Heat Mass Transf., vol. 50, no. 19–20, pp. 3749–3766, Sept. 2007. doi: https://doi.org/10.1016/j.ijheatmasstransfer.2007.02.015

[23] D. Wilcox, Turbulence Modeling for CFD, 2nd Ed., Miami, FL, USA: Amazon.com: Books, 2006.

[24] P. Crnjac, L. Škerget, J. Ravnik, and M. Hriberšek, “Implementation of the Rosseland and the P1 Radiation Models in the System of Navier-Stokes Equations with the Boundary Element Method,” Int. J. Comput. Methods Exp. Meas., vol. 5, no. 3, pp. 348–358, Abr. 2017. doi: https://doi.org/10.2495/CMEM-V5-N3-348-358

[25] Y. A. Çengel, Heat and mass transfer : a practical approach. India: McGraw-Hill Education, Pvt. Limited, 2007.

[26] G. Gordillo and H. García, Manual para el diseño y operación de hornillas paneleras. Convenio de investigación y divulgación para el mejoramiento de la industria panelera, Barbosa, Santander, Colombia: CIMPA, 1992.

[27] C. J. Greenshields, OpenFOAM User-Guide, no. 5, May. 2015.

Edxon Stiven Meneses Chacón: received his degree in Mechanical Engineer in 2014, and master’s in mechanical engineering 2018. His research interests include numerical simulation of flows involving turbulence and heat/mass transfer. Such studies are based on numerical simulations using the OpenFOAM toolkit framework. https://orcid.org/0000-0002-9413-5568

Julián Ernesto Jaramillo: received his degree in Mechanical Engineer in 1999, Master in Thermoenergetic Engineering in 2004, and Doctorate in Thermal Engineering in 2008. Professor and researcher at Universidad Industrial de Santander. https://orcid.org/0000-0001-7389-2677

Elisabet Mas de les Valls: received the Ph.D., in Nuclear Engineering in 2011 from the Polytechnic Universitat Politècnica de Catalunya (Barcelona, Spain). She is a researcher at Physics Department in the Universitat Politècnica de Catalunya. Her research interests include the study of new components for fusion technology. Such studies include magneto-convective flows, tritium transport through multiphase flows and breeding blanket design among others. https://orcid.org/0000-0003-0134-0325.