Transient Thermal Analysis in an Intermittent Ceramic Kiln with Thermal Insulation: A Theoretical Approach

Department of Mechanical Engineering, Federal University of Campina Grande, Campina Grande 58429-900, Brazil Department of Chemical Engineering, Federal University of Campina Grande, Campina Grande 58429-900, Brazil Federal Institute of Education, Science and Technology of Paraı́ba, Patos 58700-000, Brazil Department of Food Engineering, Federal University of Campina Grande, Campina Grande 58429-900, Brazil


Introduction
Research studies indicate that the first ceramic parts appeared at Dolni Vestonice (Czech archaeological site) in approximately 26000 BC, followed by Siberia in 12000 BC, China in 11000 BC, and Mesopotamia in 8000 BC, and in both Asia, the Middle East, and Europe in the late Neolithic period, between 7000 and 6000 BC [1][2][3][4]. e production of ceramic products consists of the following steps: preparation of the raw material, forming (molding), drying, and firing. In the manufacture of ceramic products, water is added to clay in order to facilitate workability in the molding step, i.e., to give plasticity to the material.
Although ceramic products have been used since ancient times, until the end of the nineteenth century, the production system of ceramic materials did not undergo major technological changes. Until this time, the production was manual, the drying was done in the sun, and the firing was carried out in field kilns, a type of intermittent kiln with low thermal efficiency [5].
In drying, an appreciable amount of thermal energy is used to evaporate the water that was added to the part during the forming period and to provide the necessary mechanical strength to reduce the chance of failure during firing process. If the drying process is not performed correctly, it can lead to material defects, loss of productivity, and increase in the energy consumption. According to the literature [6][7][8][9][10][11][12][13][14][15][16][17][18], the operating variables that influence the drying process are temperature, relative humidity and velocity of the drying air, time, shape, thickness, area/volume ratio, particle size, initial moisture content, and properties of the ceramic material.
After the drying step, the molded part is subjected to high temperatures (firing) in order to provide it rigidity and mechanical resistance. e firing process of ceramic materials in kilns must follow a preestablished firing curve. e firing curve relates the temperature of the part to the processing time, indicating the heat transfer rate at which the part must be heated or cooled and the operating temperature at each instant of time. Figure 1 shows a typical firing curve for ceramic products.
For each material, there is a critical rate of heating and cooling in such a way that if it is exceeded, product quality losses will surely occur. Optimizing the firing curve for a given material means knowing the critical heating and cooling rates in each temperature and product type.
Kiln can be defined as a structure within which it is possible to heat materials at high temperatures, i.e., it is the equipment where the drying and firing steps of the molded part occur. In the ceramics industry, a lot of equipment are made of refractory materials, which are ceramics capable of maintaining their strength and physical integrity at high temperatures (above 3000°C) and are classified into two main groups: intermittent or continuous (tunnel type) kiln [4,5,19].
Various types of fuels can be used in ceramic kilns, such as biomass (in the form of untreated wood waste and firewood), fuel oil, and natural gas [20]. Among them, natural gas stands out because it is more efficient and is considered a clean source of energy, when compared to other fuels.
Despite the great importance of an energy analysis in thermal equipment related to thermal insulators, involving the academy and the industrial sector, several works are directed toward a simple steady-state approach [21][22][23][24].
us, there is a need for more sophisticated studies, for example, those related to the transient heat transfer process.
Cavalcanti et al. [25] quantified heat losses in an intermittent kiln, concluding that convection heat loss is very low compared to radiation heat loss. e authors attribute this fact to the high surface temperature of the kiln and the low wind speeds in the surroundings (0.8 m/s).
Utlu and Hepbasli [26] performed energetic and exergetic analyses in a continuous kiln (tunnel type) for drying and firing of ceramic materials, using operational data. e energy and exergy efficiencies obtained were 39.98% and 16.41%, respectively.
Almeida et al. [27] studied numerically and experimentally the drying of hollow ceramic bricks in continuous cross-flow industrial dryer. e authors found through energy efficiency that the dryer provides a large waste of energy (over 90%) during the drying process. Low exergy efficiency was also observed, ranging from 7% to 14% depending on temperature, indicating that the drying process in this kiln is typically dissipative and therefore has high thermal energy consumption.
Gomez et al. [28] quantified the heat transfers that occur in an intermittent ceramic kiln during the heating and cooling stages. Results indicate that the greatest heat loss occurs by radiation to the sidewalls of the equipment and that a considerable amount of energy is required to heat the base, ceiling, and sidewalls of the kiln. Further, it was observed, from numerical results, that the use of 25 mm thick ceramic fiber on the sidewalls of the kiln was sufficient to provide a considerable reduction in the maximum external surface temperature and an energy gain of approximately 35% when compared with the kiln without thermal insulation.
us, in contribution to this research area, the main purpose of this work is to quantify the influence of the type and thickness of thermal insulators applied on the external sidewalls of an intermittent ceramic kiln, operating with natural gas as fuel, on the heat losses by convection and radiation, temperature distribution in the insulating material, and energy gain that the configuration provides for the system compared to the kiln without thermal insulation, using a transient approach.

Materials and Methods
is paper presents thermal analyses in an intermittent ceramic kiln, operating with natural gas as fuel, illustrated in Figure 2. e dimensions of the kiln and the experimental procedure for monitoring external ambient and kiln temperatures are reported in previous work [28].
To evaluate the influence of the type and thickness of thermal insulation on the response variables, it was considered that the temperatures on the external (T s,ext ) and internal (T int ) sidewalls of the kiln with thermal insulation during the heating step are equivalent to those obtained for the experiment performed with the kiln without thermal insulation ( Figure 3). is consideration ensures that the firing curve is the same for all analyzed cases.
With this new configuration, the temperature of the outer surface of the insulation, T s2,ext , is unknown. us, it must be determined so that it is possible to calculate the heat lost by convection and radiation to the external environment.
Four types of thermal insulators for application in the kiln were analyzed, with thickness ranging from 0.5 mm to 100 mm. To select the types of thermal insulation, factors such as market availability, temperature range, whose upper  Advances in Materials Science and Engineering limit must be greater than 300°C, and the availability of information regarding the main properties required for the calculations developed here were taken into consideration. Table 1 presents the thermal conductivity as a function of operating temperature as well as the recommended application range for each type of thermal insulation [29,30]. Table 2 indicates the density, specific heat, and emissivity of each insulating material [31][32][33][34], which do not change significantly with temperature. e general differential equation for describing onedimensional transient heat transfer is given by the following equation [35]: where ρ, c p , T, k, and S are density, specific heat, temperature, thermal conductivity, and the source term, respectively. e initial condition for thermal insulation in the sidewalls of the kiln is given by the following equation: (2) e boundary conditions on the kiln wall (x � 0) and on the outer surface of the insulation (x � L insu ) are described, respectively, by equations (3) and (4): where k insu , h r , h c , T s,ext , T s2,ext , and T amb are, respectively, the thermal conductivity of the thermal insulators, convective and radiation heat transfer coefficients, surface temperature outside the kiln, surface temperature outside the insulation, and external ambient temperature where the kiln is located.
To solve equation (1), the numerical method based on the finite volume formulation was used. e first step of the method is to discretize the domain, that is, to divide it into subdomains, also known as control volumes. e second step of the finite volume method is what distinguishes it from other numerical methods used in computational fluid dynamics: integrate the differential equation (equation (1)) under each of these subdomains, obtaining a set of algebraic equations relating the control volumes.
To understand the method, consider a one-dimensional domain that has been discretized into five control volumes, as shown in Figure 4. Taking the third control volume, its center is identified by the letter P. e east and west faces of the analyzed control volume are identified by e and w, respectively. Similarly, the centers of the east and west control volumes of this are identified by E and W. e distance between faces w and e is given by Δx, which is the length of the control volume P. e distances between the centroids P and E and between W and P are identified by δx PE and δx WP , respectively. Similarly, the distances between centroid P and face e and between face w and centroid P are identified by δx Pe and δx wP , respectively. e discretization technique was such that the west face of the first control volume coincides with boundary 1 (x � 0) and the east face of the last control volume coincides with the boundary 2 (x � L insu ), respectively, ensuring that all control volumes are whole [36]. In addition, all control volumes are of the same size (δx WP � δx PE � Δx).
Integrating equation (1) on a given one-dimensional control volume and over a time interval from t to t + Δt gives the following equation: As this analysis is applied for one transient problem, the discretization scheme can be made explicit, implicit, or fully implicit. Details of these methods can be found in Patankar   Advances in Materials Science and Engineering 3 [37] and Smith [38]. e fully implicit scheme was employed because it is unconditionally stable [35].

First Control Volume (Boundary 1).
Disregarding the source term and solving the integrals of equation (5) for the first control volume, we obtain the following equation: Taking T t+Δt P � T P and T t P � T 0 P and rearranging equation (6), the discretized equation by the finite volume method for the first control volume is given by equation (7), as follows: where k f1 and T s,ext are, respectively, the thermal conductivity and the temperature of the thermal insulation on the face in contact with the kiln wall (x � 0).

Central Control Volume.
Disregarding the source term and solving the integrals of equation (5) for the central control volumes, we have Considering linear interpolation functions for temperature and making T t+Δt P � T P and T t P � T 0 P , we can write Rearranging, we obtain the discretized equation for the central control volumes as follows:

Last Control Volume (Boundary 2).
Disregarding the source term, solving the integrals of equation (5), and considering linear interpolation functions for temperature, we have the following equation for the last control volume: e surface temperature outside the insulation, T s2,ext , can be written as a function of the heat lost by unit area of the thermal insulation (q f2 ″ ) as follows: Substituting the expression of T s2,ext obtained in equation (12) at the boundary condition of the outer surface of the insulation (equation (4)) and rearranging the terms, we obtain the heat flux lost by the thermal insulation (q f2 ″ ) as follows: Substituting the surface temperature outside the insulation (equation (12)) as a function of the heat flux lost by the thermal insulation (equation (13)) in equation (11) and rearranging the terms, we obtain the equation discretized by the finite volume method for the last control volume as follows: As can be seen, both sides of the discretized equations (equations (7), (10), and (14)) contain temperatures in the new time step, characteristic of the fully implicit discretization scheme. us, a system of algebraic equations, with one equation for each control volume, must be solved at each step of the process. For solving such systems of equations, the iterative calculation tool available in Microsoft Excel software was used. e variables k w , k e , k f1 , and k f2 are, respectively, the thermal conductivities on the west, east, west face of the first control volume, and east face of the last control volume. For the sake of simplicity, such thermal conductivity values were considered constant for all control volumes at a given time, but changing from one instant to another as a function of the average thermal insulation temperature, according to the following equation: is the average temperature of the thermal insulation at time t, T i (t) is the temperature at the center of the control volume i at time t, and N is the number of control volumes. e variables h r and h c in equation (14) are, respectively, the radiation and convection heat transfer coefficients. e radiation heat transfer coefficient is calculated according to equation (16) [22], as follows: where ε insu is the emissivity of thermal insulation, σ is the Stefan-Boltzmann constant (σ � 5.67 × 10 − 8 W/(m 2 × K 4 )), and T s2,ext and T amb are the absolute temperatures of the outer surface of the insulation ( Figure 3) and of the external surroundings, respectively. e convective heat transfer coefficient can be calculated as a function of the average Nusselt number (Nu L ), according to the following equation: where k f is the thermal conductivity of the fluid (ambient air) and L is the characteristic length, which corresponds to the height of the outer sidewall of the kiln. Modeling the kiln sidewalls as vertical flat plates, considering that the nature of the fluid flow is natural convection, and admitting laminar flow (10 4 ≤ Ra L ≤ 10 9 ), the correlation proposed by Churchill and Chu [39] was used to calculate the average Nusselt number, as follows: where Pr is the Prandtl number and Ra L is the Rayleigh number, calculated from equation (19), as follows: where ] is the kinematic viscosity, α is the thermal diffusivity, g is the gravity acceleration, β � 1/T f is the volumetric expansion coefficient, T s2,ext and T amb are the outer surface temperature of the thermal insulation and external ambient temperature, respectively, and L is the height of the outer sidewall of the kiln. e values of the parameters ], α, Pr, and β are tabulated according to the fluid type and the film temperature (T f ), which corresponds to the average temperature between the external ambient temperature (T amb ) and the outer surface temperature of the thermal insulation (T s2,ext ), and their values are calculated for each time of the process.
To guarantee that the numerical results were independent of the number of control volumes, a mesh convergence analysis was developed using the methodology proposed by Celik et al. [40], which is based on Richardson's extrapolation [41,42]. Using this methodology, it is possible to determine the ideal mesh by calculating the grid convergence index (GCI). In addition, Richardson extrapolation allows to estimate the exact solution from the results of existing meshes.
is methodology has been used in several works related to the most diverse areas of computational fluid dynamics [43][44][45][46][47][48]. e entire procedure for mesh convergence analysis is described in detail in Gomez et al. [28].
In addition to the mesh convergence analysis, a time step independence study based on absolute error was employed. e response variable (ϕ) for the two previously mentioned studies was the energy gain by the thermal insulation during the kiln heating process (Q in insu ), calculated using the composite trapezoidal rule, as follows: where

Advances in Materials Science and Engineering
where a is the time instant that the kiln was turned on (a � 0 s), b is the time instant that the kiln was turned off (b � 840 min, based in experiments), and n is the number of time steps in the heating step. q in insu (a), q in insu (b), and q in insu (i) are the rates of energy entering the thermal insulation at time a, b, and i, respectively.
To quantify the influence of type and thickness of thermal insulation on the kiln efficiency, it was proposed the energy gain variable, which can be calculated according to the following equation: where E supplied and E supplied insu are the energies provided during the heating process for the uninsulated and thermally insulated kiln, respectively, and are calculated according to equations (23) and (24), as follows: where _ E supplied and _ E E supplied insu are the heat rates supplied to the kiln without and with thermal insulation, respectively, and are calculated according to equations (25) and (26), as follows: where _ E st is the rate at which energy is stored in the kiln, q lost (base/ceiling) is the heat lost by the base/ceiling assembly of the kiln, and q lost sidewalls is the heat lost by the sidewalls of the kiln without thermal insulation. e calculation methodology for each heat transfer parameter in the kiln without thermal insulation is presented in detail in Gomez et al. [28]. q lost insu is the heat lost by thermal insulation (equation (27)) and _ E st insu is the rate at which energy is stored in the thermal insulation, calculated according to equation (28).
where A is the sum of the kiln sidewall areas; ρ insu , L insu , and c p insu are, respectively, the density, thickness, and specific heat of the thermal insulators; and ΔT average insu is the average temperature variation within the thermal insulation in the time interval ∆t.
It is important to emphasize that, given the considerations that have been made (Figure 3), the curves of q lost(base/ceiling) and _ E st as a function of time are common for the kiln with and without thermal insulation. What changes from case to case is the heat lost by the sidewalls of the kiln, influenced by the type and thickness of the thermal insulation. Figure 5 shows the average temperatures of the outer and inner surfaces of the kiln without insulation, as well as the ambient air temperature in the surroundings of the kiln during the heating step, obtained by experiments. It is observed that the curve for the external surface temperature of the kiln, T s,ext , is one of the boundary conditions used in the analysis of the kiln with thermal insulation, as seen in equation (3). e curve for the ambient air temperature outside the kiln (T amb ) is one of the input parameters of the other boundary condition of the physical problem, as seen in equation (4).

Results and Discussion
As previously mentioned, a mesh convergence study was performed representing the thermal insulation used in the kiln sidewalls. For each thermal insulation thickness analyzed, three types of meshes were created, named M 1 , M 2 , and M 3 , with 20, 10, and 5 control volumes, respectively. Table 3 presents the representative mesh size as a function of the thermal insulation thickness, where l 1 , l 2 , and l 3 are the representative mesh sizes of M 1 , M 2 , and M 3 , respectively. In order to satisfy the recommendations proposed by Celik et al. [40], the number of control volumes of each mesh was defined such that the refinement factors between the meshes were equal, i.e., r 21 � r 32 � 2. It is worth noting that all analyses were performed using the same type of thermal insulation (ceramic fiber) and the same time step (Δt � 0.1 min). e response variable used to do the mesh convergence analysis was the heat conduction energy gain by the thermal insulation during the heating stage (Q in insu ), according to equation (20). e values of this variable as a function of the mesh and thermal insulation thickness are reported in Table 4. Table 4 presents the results of the mesh convergence study parameters as a function of thermal insulation thickness (L insu ).
e values of C in the table indicate monotonic convergence for cases where L insu ≥ 7.5 mm, since their values are in the range between 0 and 1, ensuring the applicability of the Richardson extrapolation method to the variable of interest in the given range [41]. For L ≤ 5.0 mm, it is observed that the convergence is oscillatory, since − 1 < C < 0. Celik et al. [40] states that the oscillatory convergence should not be considered as an unsatisfactory result, because if ε 21 � ϕ 2 − ϕ 1 or ε 32 � ϕ 3 − ϕ 2 is very close to zero, it may indicate that the exact solution has already been reached, as it seems to be the case. It is important to emphasize that the step-by-step procedure to calculate the C parameter and the mesh convergence index (GCI) can be found in the previous work [28].
It is observed that there is a reduction in the convergence condition for all analyzed cases, since GCI 21 < GCI 32 , indicating that the dependence was reduced with the mesh refining. From the proximity of the values of GCI 32 and (r 21 ) p GCI 21 , it is possible to state that the asymptotic interval was reached and that the extrapolated solution is close to the exact solution for this variable, for all analyzed cases [49]. For all thickness values, there is also a good proximity between the extrapolated solution and the solutions for M 1 and M 2 . us, considering that the more refined the mesh, the greater the total simulation time, the 10 control volumes mesh (M 2 ) was chosen for further analysis.
For all analyzed cases, the energy balance is not satisfied for the 5-volume control meshes. e 10 and 20 control volumes have values below 10 − 9 for the global energy balance in thermal insulation during all time instants of the heating process, considering a convergence criterion of 10 − 7 . e next step was to perform a time step independence analysis using the same type of thermal insulation (ceramic fiber), the same mesh (M 2 ), and the same response variable (Q in insu ). For this analysis, three distinct time steps were used, as presented in Table 5. It is observed, for all analyzed cases, that the absolute error of the response variable between the time steps of 0.1 min and 1.0 min (ε 12 � ϕ 1 − ϕ 2 ) is considerably smaller when compared to the absolute error between the time steps of 1.0 min and 15.0 min (ε 23 � ϕ 2 − ϕ 3 ). us, considering that the smaller the time step used, the greater the total simulation time, a time step of Δt � 1 min was adopted for further analysis. Figure 6 illustrates the average temperature in thermal insulation during the heating process for several ceramic fiber thicknesses. ese values are calculated from the arithmetic mean of the temperatures obtained in the center of each of the 10 control volumes at each time instant.
Such results are essential for the analysis, since the thermal conductivity of the insulating material for each control volume at a given time is iteratively calculated as a function of the average temperature in the insulation at its respective time instant (equation (15)). From the analysis of this figure, it was verified that regardless of the thickness of the thermal insulation, the average temperature tends to increase along the heating process, as expected. As seen in Table 1, the thermal conductivity of the insulating materials increases with the increase of the average temperature in it; thus, the increase in the insulating material thickness promotes a decrease in the average temperature of the thermal insulation and, consequently, a reduction in its thermal conductivity. Figure 7 illustrates the external surface temperature of the thermal insulation during the heating process for several ceramic fiber thicknesses.
From the analysis of this figure, it is observed that the external surface temperature of the thermal insulation decreases with an increase in the material thickness, proving that the greater the thickness of the insulating material, the higher the thermal resistance to heat transfer. us, a higher thickness of thermal insulation implies a smaller difference between the external surface temperature and the ambient temperature, thus promoting a reduction in convection and radiation heat losses through the kiln sidewalls, as will be confirmed below. It is important to state that the 0 mm thickness curve in Figure 7 represents the temperature of the thermal insulation on the face in contact with the kiln wall, as seen in Figure 3 and equation (3), therefore being the first boundary condition (x � 0) of the physical problem for the thermally insulated kiln. Figure 8 shows the convective heat transfer coefficient during the heating process for several ceramic fiber thicknesses. e thermophysical properties, kinematic viscosity, thermal diffusivity, volumetric expansion coefficient, and

Advances in Materials Science and Engineering
Prandtl number, tend to increase the value of the convective heat transfer coefficient as the thermal insulation thickness is increased. is happens because the increase in the thermal insulation thickness promotes a reduction in the surface temperature outside the insulation (Figure 7) and, consequently, a reduction in the film temperature. In contrast, the increase in the thermal insulation thickness promotes a decrease in the thermal conductivity of the fluid (k f ), because of the reduction in the film temperature [22], and in the difference between the surface temperature outside the thermal insulation and the ambient temperature, as previously seen, both contributing, more significantly, to a reduction in the convective heat transfer coefficient. Figure 9 shows the radiation heat transfer coefficients during the heating process for several ceramic fiber thicknesses. e radiation heat transfer coefficient strongly depends on the temperatures of the outer surface of the thermal insulation and of the external surroundings, as seen in equation (16), making its value more influenced by the thermal insulation thickness as compared to the convection heat transfer coefficient. Figures 10 and 11 show the heat fluxes lost by convection and radiation, respectively, during the heating process, for several ceramic fiber thicknesses.
As expected, the greater the thermal insulation thickness on the kiln sidewalls, the lower the heat lost by convection   and radiation. It is also noted that for a given thermal insulation thickness, the heat lost by radiation from the sidewalls is always greater than the heat lost by convection. is is because the radiation heat transfer coefficient is greater than the convection heat transfer coefficient throughout the kiln heating process. Figure 12 illustrates the rate at which energy is stored in the thermal insulation during the heating process for several ceramic fiber thicknesses.
e analysis of this figure shows that the greater the insulation thickness is, the higher the rate at which energy is stored in it. Among the four types of thermal insulation analyzed, ceramic fiber has the highest accumulated energy rates. is is because both the density and the specific heat of the ceramic fiber are higher than those presented by the other insulating materials analyzed (Table 2), making it have a greater capacity to absorb heat. Figure 13 shows the heat rate that must be supplied to the kiln during the heating process for several ceramic fiber thicknesses.
From the analysis of this figure, it can be seen that up to approximately 100 minutes, the heat transfer rates provided to the kiln are very close. From this moment on, the heat flux supplied to the kiln is strongly influenced by the thermal insulation thickness. e greater the thickness of the thermal insulation is, the lower the energy rates that must be supplied to the kiln at each instant of time in order to maintain the same firing curve during the heating step. is can be explained by the fact that the increase in the thermal insulation thickness increases the thermal resistance by conduction, which reduces the surface external temperature of the thermal insulation (Figure 7), providing an increase in the convection and radiation thermal resistances. ese  phenomena reduce heat loss and the energy rate that must be supplied to the kiln. Figure 14 shows the accumulated energy supplied to the kiln during the heating process for several ceramic fiber thicknesses.
Such curves are obtained by integrating in time the heat rate curves (Figure 13) for each type and thickness of thermal insulation throughout the heating process. It is evident that the greater the thickness of the thermal insulation is, the lower the thermal energy that must be supplied to the kiln during the heating step. Figure 15 shows the influence of type and thickness of thermal insulation on the total energy supplied during the kiln heating step. From this result and using equation (22), it was possible to calculate the energy gain as a function of the type and thickness of the thermal insulation, as shown in Figure 16.
A number of observations can be made from the analysis of Figures 15 and 16 as follows. e greater the thermal insulation thickness is, the lower the energy that must be supplied during the kiln heating process and, consequently, the greater the energy gain. It is important to emphasize that from a certain thickness value (50 mm approximately), the energy gain increases insignificantly with the increase in the thermal insulation thickness.
Fiberglass is the type of thermal insulation that provides the greatest energy gain for any thickness analyzed, followed by rockwool, calcium silicate, and finally ceramic fiber.
While 5.0 mm of fiberglass is sufficient to provide an energy gain of approximately 29.9%, the same thickness for rockwool, calcium silicate, and ceramic fiber provides, respectively, energy gains of 28.0%, 26.5%, and 24.3%.
For low insulation thicknesses, the greater the influence of the insulation type on the energy gain, while for high thicknesses the energy gain is little influenced by the type of thermal insulation. For a thickness of 2.5 mm, the range of energy gain is from 17.9% (ceramic fiber) to 24.1% (glass fiber), i.e., a variation of 6.2%. For a thickness of 100 mm, the energy gain ranges from 38.3% (ceramic fiber) to 39.7%, i.e., a 1.4% variation.
e results presented here indicate that the lower the thermal conductivity of the material, the greater the energy gain, as expected. e use of a thermal insulator with lower thermal conductivity implies less heat loss and therefore less energy must be supplied to the kiln, resulting in a higher energy gain. Figure 17 shows the influence of thermal insulation type and thickness on the maximum external surface temperature obtained during the kiln heating process.
A number of observations can be made from the analysis of Figure 17, as follows. e greater the thickness of the thermal insulation, the lower the maximum external surface temperature of the kiln.
Fiberglass is the insulating material which provides the greatest reduction in the maximum external surface temperature for any thickness analyzed, followed by rockwool, calcium silicate, and finally ceramic fiber. While 2.5 mm of fiberglass is sufficient to reduce the maximum external surface temperature of the kiln from 249.3°C (without thermal insulation) to 148.1°C, the same 2.5 mm of rockwool, calcium silicate, and ceramic fiber reduces the maximum external surface temperature to 150.9°C, 158.2°C, and 176.8°C, respectively.
For thicknesses above 5.0 mm, the higher the thickness of the thermal insulation, the lower the influence of the   thermal insulation type on the maximum external surface temperature. For a thickness of 5.0 mm, the range of the maximum external surface temperature is 115.7°C (fiberglass) to 145.2°C (ceramic fiber), i.e., a range of 29.5°C. For a thickness of 100 mm, the maximum external surface temperature range is 42.8°C (fiberglass) to 49.3°C (ceramic fiber), i.e., a range of 6.5°C. e results indicate that the lower the thermal conductivity of the material, the lower the maximum external surface temperature.
In regard to fiberglass, it is a fiber-reinforced polymer using glass fiber or glass reinforced plastic [50]. It is made of plied, spun, and texturized continuous glass filaments. e fiberglass has many advantages, when compared to other types of thermal insulation, such as low thermal conductivity, low density, very high strength, high sound absorption, flame retardant (no burn), exceptional weathering resistance (low moisture absorption), and low cost [50][51][52]. us, fiberglass has been applied as thermal and acoustic insulation at different thermal equipment, pipelines, and building constructions, reducing heat transfer. It is also commonly used to prevent condensation of water vapor on cold surfaces, such as air conditioning ducts (ductwork) and cold-water pipes. Furthermore, fiberglass is not generally considered to be dangerous, but if not properly sealed off, it can get into air vents and circulate through the environment, irritating the skin and respiratory system of humans.

Conclusions
Considering what has been presented, it is important to use thermal insulation in processes where it is desired to reduce heat losses to the external environment and, consequently, to promote an increase in its thermal efficiency.
To ensure the reliability of the numerical results obtained, mesh convergence and time step analyses were performed representing thermal insulation, based on Richardson extrapolation, grid convergence index (GCI), and absolute error calculations.
From the obtained results and discussions made in this work, it can be concluded that (a) e proposed methodology can be easily applied in similar processes, having knowledge of the temperature and air velocity conditions in the surroundings of the kiln, as well as the external and internal surface temperatures of the kiln. (b) e mathematical procedure developed in this paper presented a low computational effort when compared to commercial CFD solvers, such as STAR CD ® , Ansys Fluent ® , and Ansys CFX ® .
(c) Results indicated that the higher the thickness of the thermal insulation, the lower the maximum external surface temperature and the greater the energy gain when compared to the kiln without thermal insulation. While 1 mm fiberglass thickness provides a 57.88°C reduction in the maximum external surface temperature and an energy gain of 15.37%, a 100 mm fiberglass thickness provides a reduction of 206.59°C in the maximum external surface temperature and an energy gain of 39.67% when compared to the kiln without thermal insulation. (d) Of the four types of thermal insulation analyzed, fiberglass is the insulating material that provides the greatest energy gain and, consequently, promotes the greatest economic and environmental gains. (e) Fiberglass, due to excellent physical characteristics, such as low thermal conductivity and density, is the thermal insulation that provides the greatest reduction in maximum external surface temperature, contributing to the greatest reduction in thermal discomfort and the risk of work accidents when in operation.
Finally, the results presented here can be used as a decision-making tool in choosing the type of thermal insulation and its thickness that are capable of providing the desired energy gain for the kiln, with a better benefit/cost ratio. It is important to emphasize that other factors must be considered when choosing the type of insulation, such as shelf life, maintenance cost, and effects on human health.
Abbreviations a E , a P , a 0 P , a W , S u , S P : Coefficients (-) A: Area (m 2 ) c p : Specific heat (J kg − 1 K − 1 ) c pinsu : Specific heat of the thermal insulation (J kg − 1 K − 1 ) C: Mesh convergence parameter (-) _ E st : Rate at which energy is stored in the kiln (W) _ E st insu : Rate at which energy is stored in the thermal insulation (W) _ E supplied : Heat rate supplied to the kiln without thermal insulation (W) E supplied : Energy supplied during the heating process for the kiln without thermal insulation (J) _ E insu supplied : Heat rate supplied to the kiln with thermal insulation (W) E supplied insu : Energy supplied during the heating process for the kiln with thermal insulation (J) g : Gravity acceleration (m s − 2 ) h c : Convection heat transfer coefficient (W m − 2 K − 1 ) h r : Radiation heat transfer coefficient (W m − 2 K − 1 ) GCI: Grid convergence index (-) k e : ermal conductivity on the east face of the control volume (W m − 1 K − 1 ) k f : ermal conductivity of the fluid (ambient air) (W m − 1 K − 1 ) k f1 : ermal conductivity on the west face of the first control volume (W m − 1 K − 1 ) k f2 : ermal conductivity on the east face of the first control volume (W m − 1 K − 1 ) k insu : ermal conductivity of thermal insulation (W m − 1 K − 1 ) k w : ermal conductivity on the west face of the control volume (W m − 1 K − 1 ) l: Representative mesh size (mm) L: Characteristic length (height of the outer sidewall of the kiln) (m) L insu : ermal insulation thickness (mm) N: Number of control volumes (-) Nu L : Average Nusselt number (-) Pr: Prandtl number (-) q lost(base/ceiling) : Heat lost by the base/ceiling assembly of the kiln (W) q lost sidewalls : Heat lost by the sidewalls of the kiln without thermal insulation (W) q in insu : Rate of energy entering the thermal insulation (W) q conv insu : Heat flux lost by convection (W) q rad insu : Heat flux lost by radiation (W) Q in insu : Energy gain by the thermal insulation during the kiln heating stage (J) Ra L : Rayleigh number (-) S: Source term (W m − 3 ) t: Time (s) T: Temperature (°C) T average insu : Average temperature of the thermal insulation (°C) T int : Surface temperature inside the kiln (°C) T s,ext : Surface temperature outside the kiln (°C) T s2,ext : Surface temperature outside the thermal insulation (°C) T amb � T ∞ : Ambient air temperature outside the kiln (°C) T f : Film temperature (°C) T E : Temperature at the control volume E at time t + ∆t (°C) T P : Temperature at the control volume P at time t + ∆t (°C) T 0 P : Temperature at the control volume P at time t (°C) T W : Temperature at the control volume W at time t + ∆t (mm) δx PE : Distance between the centroids P and E (mm) δx WP : Distance between the centroids W and P (mm) δx Pe : Distance between the centroid P and face e (mm) δx wP : Distance between the face w centroid P (mm) Δx: Length of the control volume P (mm) α: ermal diffusivity (m 2 s − 1 ) β: Volumetric expansion coefficient (K − 1 ) ε: Emissivity (-) υ: Kinematic viscosity (m 2 s − 1 ) ρ: Density (kg m − 3 ) σ: Stefan-Boltzmann constant (W m − 2 K − 4 ) ϕ: Response variable (MJ) ϕ 21 ext : Extrapolated solution for the response variable used in the mesh convergence study (MJ).

Data Availability
e data used to support the findings of this study are included within the article.