Thermal properties determination of a cylindrical product during its cooling: two-dimensional numerical model and uncertainty

ABSTRACT In this article, a two-dimensional model is proposed to determine thermal diffusivity and convective heat transfer coefficient, providing the average values, their uncertainties and the covariance matrix referring to these parameters for a product with cylindrical geometry during its cooling. The proposed model used a two-dimensional numerical solution of the diffusion equation for the direct problem; an experimental dataset referring to the cooling kinetics in the center of the product and an optimizer program based on the Levenberg-Marquardt algorithm for the inverse problem. The model was used with success to determine the thermal properties, their uncertainties and the covariance matrix of a cucumber during its cooling. The obtained results allowed establishing a confidence band that made it possible to graphically evaluate the precision of a new simulation for a cucumber with different dimensions.


Introduction
In order to increase the post-harvest life of agricultural products, their processing is usually carried out, often involving the simultaneous heat and mass transfer, as in the case of drying. [1][2][3] Other types of processing involve simple heat transfer, such as pasteurization [4][5][6] , cooling [7][8][9][10][11] , and also freezing. [12,13] To simulate all these processes, and consequently to obtain information such as their duration, as well as the energy required, among others, it is necessary to know the initial and final temperatures of the product, its geometry, and its dimensions, as well as its thermophysical properties. To determine these properties during a transient state, various optimization techniques can be used such as the concept of lag factor [7] or the similar "slope method". [10] This technique uses a single term representing the analytical solution of the partial differential equation which describes the physical phenomenon and it disregards, by graphical inspection with the use of logarithm, the first experimental points, in such a way that this single term is sufficient to describe well the phenomenon. However, if many experimental points are disregarded, statistical information is lost. Thus, Silva et al. [8] proposed an algorithm for the optimal determination of the number of points to be disregarded. To this technique of determination of parameters, the authors gave the name of "Optimal Removal of Experimental Points -OREP", and such technique makes it possible to generate good results for the thermal properties during the transient state of a process.
Some methods (so-called "robust"), which determine parameters of partial differential equations through successive trials and an experimental dataset, are also available in the literature. [3,6] However, for the specific case of heat transfer, in many works available in the literature the values of the thermal properties are determined without their uncertainties being known, as well as the covariance between such values. Thus, the use of the values of such properties to simulate pasteurization, cooling, and freezing are compromised due to the lack of information on the precision of the results obtained.
Le Niliot, and Lefevre [14] and also Mariani et al. [12] , in order to determine the uncertainties of thermal properties, proposed to perform 100 numerical simulations where the measurements are disrupted with 100 different Gaussian distributions (with zero mean and standard deviation obtained from the simulation using the parameters originally determined). Thus, the uncertainties were obtained by performing the statistical treatment of the 100 values for each parameter. Da Silva and Silva [9] have also used this method to determine the covariance matrix related to the thermal diffusivity and the convective heat transfer coefficient during cucumber cooling. However, this method, despite being able to determine the uncertainties of parameters and the covariance matrix, is very slow due to the number of simulations required. Thus, this method has usually been employed for one-dimensional geometries, whose simulations are generally faster than those for two-and three-dimensional ones.
Other techniques to determine parameters of differential equation (and their uncertainties) using experimental dataset have also been found in the literature in recent years. [15][16][17] However, the necessary codes to do it (such as OptiPa, SBtoolbox2, and AMIGO2) are usually made available (or developed) for MATLAB environment. In the year 2017, the first two authors of the present article developed the program LS Optimizer (Version 6.2, 2018, http://zeus.df.ufcg.edu.br/labfit/LS.htm), to be installed directly on Windows platform. This free software, based on the Levenberg-Marquardt algorithm [18][19][20] , is very fast to determine parameters of a differential equation using experimental data. Thus, this software makes it possible to use two-and three-dimensional models in order to determine thermal properties during a transient state, with no need for a specific computational environment. In this context, the objectives of this article are defined below.
The main objective of this article was to propose a two-dimensional numerical model to determine thermal properties (their uncertainties and the covariance matrix) of products with cylindrical geometry, using an experimental dataset obtained during their cooling. Thus, the obtained values can be used to predict new cooling curves for the same product, but with different dimensions. As an application, the model was used to determine the thermal diffusivity and the convective heat transfer coefficient (their uncertainties and covariance matrix) for a cucumber during its cooling process.

Heat conduction in cylindrical geometry
If a heat conduction process in a cylindrical domain (with radius R and height L) is symmetric with respect to the y-axis, as shown in Figure 1a, the diffusion equation to describe the unsteady state process is given by: where ρ (kg m −3 ) is the density, c p (J kg −1 K −1 ) is the specific heat, T (K) is the temperature, k (W m −1 K −1 ) is the thermal conductivity, t (s) is the time and r(m) and y(m) are the coordinates of position. For agricultural products, a cooling process usually involves a temperature variation of about 20°C or less. Thus, the following assumptions can be accepted to describe cooling of cylindrical agricultural products: (1) the thermophysical parameters can be considered as constant during the process; (2) for some of these products, it can be assumed that the cylindrical medium is homogeneous and isotropic; (3) for the proposed model, it was considered that heat transfer can be described only by conduction. Thus, the thermal diffusivity α (m 2 s −1 ), defined below, should be considered as "apparent": Since the thermophysical properties were considered as constant, Equation (1) can be divided by ρc p , and rewritten as:

Direct problem: solver
A solver for Equation (3) was obtained using the Finite Volume Method with a fully implicit formulation. [21] Several authors have shown that the diffusion equation for the cylinder has unique analytical solution. [8,11] However, for the numerical solution proposed, the fully implicit formulation was chosen because it is unconditionally stable. [9] A numerical solution was chosen because such a solution can be obtained even if the thermal properties and the dimensions of the domain are variable along the process. On the other hand, considering a symmetric heat conduction with respect to the y-axis, only the piece shown in Figure 1b can be considered to solve the diffusion equation. Thus, a grid can be created in the rectangle presented in Figure 1c. As a result, a two-dimensional uniform grid can be obtained, as shown in Figure 2a. In addition, Figure 2b shows a rectangular  element of the grid, whose sides are Δr and Δy. A revolution of this rectangular element in Δθ radians around the y-axis shows the control volume given in Figure 2c: ΔV ¼ Δθ r P Δr Δy, where r P defines the radial position of the nodal point P of the rectangular element. By observing Figure 2a, it is possible to conclude that, for the two-dimensional domain, there are nine types of control volumes: internal (no contact with external medium), north, north-east, south, south-east, east, north-west, south-west, and west. Due to the radial symmetry, it is observed that the last three types of control volumes have flux zero at the west boundary. Integrating Equation (3) in space (Δθ r P Δr Δy) and time (Δt) for a given control volume, after the simplifications, yields: @T @r e À r w α w @T @r w Δy þ α n @T @y n À α s @T @y where the superscript 0 means "former time" t, and its absence means "current time", t+Δt. The subscripts "e", "w", "n", "s" are, respectively, east, west, north, and south interfaces, while P is the nodal point of the control volume.

Internal control volume
A control volume is defined as "internal" if it has neighbors to the north, south, east, and west. As an example, an internal control volume and its neighbors to the north, south, east, and west are shown in Figure 3a. For this type of control volume, substituting the derivatives of Equation 4 by numerical approximations of first-order yields: where: A n ¼ r P α n Δr Δy (9) East boundary A control volume P at the east boundary and its neighbors to the north (N), south (S), and west (W) are shown in Figure 3b. In this figure, the symbols denoted by T 1;e (K) and h e (ms −1 ) represent the equilibrium temperature and the convective heat transfer coefficient at the east boundary, respectively. Thus, the following result for the discretization is obtained: The coefficients A w , A n , and A s of Equation 12 are given by Equations 8-10, respectively. The coefficients A P and B are given by: and As an additional information, this numerical solution has been previously validated using a twodimensional analytical solution presented in Ref. [11] On the other hand, the last results were obtained by imposing equal diffusive and convective fluxes at the east interface of the control volume with nodal point P (see Figure 3b): À α e @T=@r e ¼ h e ðT e À T 1; e Þ. This equation defines the boundary condition of the third kind for the east side (subscript "e"). In this last equation, T e (K) is the temperature at the east boundary, while T 1; e is the cooling temperature at the east side of the medium. On the other hand, the heat transfer coefficient h He (W m −2 K −1 ) and h e are related as follows: The same discretization procedure mentioned above can be used to discretize the other seven types of control volumes. Thus, a system of equations for each time step is obtained with the Equations 5,12 and other similar equations, and this system can be solved, for instance, by the Gauss-Seidel method [22] , with a convergence tolerance of 1 ×10 −8 . Specifically in this article, because the thermal properties are considered constant, for each control volume the following equality is observed: Also, at the boundaries of the rectangular grid ( Figure 2a) the following equality can be written: On the other hand, on the whole west side, due to the radial symmetry, the heat flux is zero. The solver for the simulation of heat transfer at a specified point within the finite cylinder was created in FORTRAN. This solver was created in Compaq Visual Fortran Professional Studio, Edition V. 6.6.0, using a programming language option called QuickWin Application.

Inverse problem: parameter determination
Thermal diffusivity α and convective heat transfer coefficient h were determined using LS Optimizer Software, developed by the first and second authors of this article. The program is freely available at http://zeus.df.ufcg.edu.br/labfit/LS.htm, and it is intended for the solution of inverse problems, determining parameters of differential equations and other types of complex functions for a specified experimental dataset. The Levenberg-Marquardt algorithm [18][19][20] is used by LS Optimizer, which executes the solver provided by the user in order to determine the parameters that allow a simulation as close as possible to the experimental dataset. The software's Help Menu explains how to add a few lines in the solver's original code, allowing the reading of information provided by LS, and vice-versa, through text files.
According to the information available in the LS Optimizer manual, as long as the solver for the direct problem and the experimental dataset are provided, this optimization software executes the solver in order to obtain the necessary information for the determination of parameters, using the statistical indicator chi-square as the objective function. As mentioned by Da Silva et al. [11] , through the experimental (T exp i ) and simulated (T sim i ) values, the optimizer program iteratively determines corrections for h and α, and calculates the chi-square referring to each simulation through the expression =ðN À nÞ q , in which (N-n) is the number of degrees of freedom. Hence, at the end of the iterative process, one can impose σ T i ¼ σ and recalculate the covariance matrix, which provides the uncertainties of the parameters and the correlation between them. Thus, starting from initial values α 0 and h 0 , the optimizer corrects these parameters iteratively up to convergence and, in the end, LS Optimizer delivers the optimal values for α and h, their uncertainties, the covariance matrix, and graphs.
In this article, the initial values for the parameters to be determined were established as α 0 = 1 x 10 −7 m 2 s −1 and h 0 = 1 x 10 −6 m s −1 . If the LS user considers that the initial values are not so reasonable, the corrections calculated by the optimizer can produce the divergence of the iterative process that obtains the parameters. Due to that, the optimizer has a resource called "split first corrections into…". In this optimization process, using this resource, the first corrections of the parameters were divided into 40 parts. Therefore, only 1/40 of the values calculated were used to correct the parameters in the first corrections.

Experimental dataset
The model proposed in this article to determine thermal properties of a cylindrical product involves: (1) A solver for the direct problem, presented in Section 2.2; (2) An optimizer program for the inverse problem, presented in Section 2.3; and (3) An experimental dataset referring to the cooling kinetics at a specified point within the cylindrical domain.
Silva et al. [8] have proposed a one-dimensional model in which the thermal properties are determined by the fitting of only the first term of the series representing the analytical solution of the onedimensional diffusion equation. Due to this proposal, the first experimental points must be disregarded to perform the curve fitting. Then, the authors have proposed an algorithm to determine the number of the optimal removal of experimental points, in order to maintain as much statistical information as possible. This model was applied to the cooling of cucumbers under natural convection and, in the present article, for comparison purposes, the same experimental dataset will be used. The cucumber used in the experiment has a radius R = 0.019 m and a height L = 0.160 m. During the experiment, the mass loss was considered to be negligible. The moisture content of the product was X = 0.96 (w.b.), and its initial temperature was T 0 ¼ 22:0 o C. The temperature of the cooling air was T 1 ¼ 4:0 o C, and its velocity was kept at 2 m s −1 , with a relative humidity of 80%. The cooling kinetics was determined through a thermocouple placed in the center of the product (r = 0 and y = 0), and the temperature was measured at regular intervals of time, totaling 37 experimental points. In order to determine the thermal properties, the dataset was written in the dimensionless form:

Results and discussion
In order to solve the direct problem, a 40 × 100 grid was used and the cooling duration was divided into 1000 time steps. For these values, it was possible to observe that the first six significant figures of each temperature obtained are independent of new refinements for the grid and time.

Results
After convergence, LS Optimizer provided the values of the parameters, and also the factor to multiply the uncertainties so that the confidence interval is 95.4% (in the present case, it was obtained: factor = 2.04). Thus, the final results for the thermal properties and their uncertainties are given in Table 1.
It is interesting to observe that the thermal conductivity k was determined from Equation 2, and the heat transfer coefficient h H was calculated using Equation 15. For that, the specific heat was estimated from the expression [23,24] : c p ¼ 1:381 þ 2:930X ðkJ kg À1 K À1 Þ. The cucumber density value used in the present article was determined by Fasina and Fleming [25] as: ρ = 959 kg m −3 . It is interesting to observe that the value obtained for k agrees with the value obtained by the empirical equation provided in Ref. [23] for the thermal conductivity (k = 0.148 +0.493X): k = 0.62 W m −1 K −1 . The discrepancy between the two values is only 3%.
Another result provided by LS Optimizer is the covariance matrix, given by:, cov ¼ 3:672 Â 10 À17 À 7:021 Â 10 À16 À7:021 Â 10 À16 1:491 Â 10 À14 ! indicating a correlation coefficient between α and h of −0.9489 and, even with this high negative correlation, the iterative process has converged. With the results obtained for α and h, the cooling kinetics of the cucumber is shown in Figure 4. The statistical indicators chi-square and determination coefficient for the result presented in Figure 4 are χ 2 = 2.6981 ×10 −3 and R 2 = 0.9991, respectively. The error distribution referring to the cooling kinetics is presented in Figure 5, and the residuals were calculated through the differences between experimental and simulated dimensionless temperatures, T Ã exp i À T Ãsim i for all experimental points i. In Figure 5, the continuous line represents the average error, and its value was 9.48 ×10 −4 , in good agreement with the value zero, which is theoretically expected.

Simulations for a cucumber with other dimensions
Obviously, through the results obtained in this article, it is possible to simulate cooling curves for cucumbers under similar conditions, but with other dimensions, with no need for new experiments. As an example, for a cucumber with a radius R = 0.026 m and a height L = 0.220 m, the cooling kinetics at the central point, shown in Figure 6, can be obtained by simulation. On the other hand, Table 1. Thermal properties (thermal diffusivity, thermal conductivity, convective heat transfer coefficient, heat transfer coefficient) determined by the proposed model and by Silva et al. [8] Thermal property Proposed model Silva et al. [8] Figure 6a shows that the confidence band is very narrow, indicating that the results obtained for α and h have a good precision to simulate the transient state of a cooling process. Figure 6b shows that, after 4320 s, the most probable value of the dimensionless temperature at the center of the second cucumber is 0.220 (about 8.0°C), while for the original cucumber a value of 0.102 (about 5.8°C) was experimentally found, as shown by Figure 4.

Discussion
In the literature, in order to simulate a transient state, many researchers determine one thermal property through empirical correlations. [10,15] Due to that only one thermal property is determined using some optimization algorithm and an experimental dataset. In the model proposed in this article, two thermal properties (α and h) were simultaneously determined by optimization, which allowed determining the covariance between these properties. Thus, the proposed model allowed establishing a confidence band for cooling simulations of the product, using the obtained results. On the other hand, as it was observed, the proposed model included a numerical solution of the twodimensional diffusion equation for the cylinder. Consequently, this model is useful not only to describe cylindrical product cooling, which involves a temperature variation of about 20°C or less, and therefore the thermal properties can be considered constant [10,11,[26][27][28] ; but also to determine thermal properties during, for instance, the pasteurization of a product, which involves a large temperature variation and, therefore, such properties can be considered variables. [6] Due to the use of LS Optimizer Software, the main statistical indicators for each obtained parameter can be provided. In addition, according to Silva et al. [3] , citing Taylor [29] , the Student's t-test is a statistical indicator that allows determining the probability Pð V=σ V Þ that the average value V of a parameter given in the form V ¼ V Ç σ V is zero, despite its value determined by optimization. In order to obtain this information for α and h using this test, the values of Pð V=σ V Þ were calculated, and both results were equal to zero.
Silva et al. [8] , using a one-dimensional analytical model and a technique called OREP, has obtained α = 1.47 ×10 −7 m 2 s −1 and h = 6.39 ×10 −6 ms −1 for the same dataset used in this article. These values are compatible with the two-dimensional numerical model proposed in the present article, which delivered the results α = (1.48 ± 0.12)×10 −7 m 2 s −1 and h = (6.35 ± 0.25)×10 −6 m s −1 , with 95.4% confidence level. Although the two models have generated compatible results, the proposed model also provided the covariance matrix, which made it possible to determine not only the uncertainties of the thermal properties that were calculated but also the correlation between them. As mentioned before, this information is important because it allows determining a confidence band for the cooling kinetics of the product, as shown in Figure 6a. As an example, for t = 1441 s, it is possible to estimate the dimensionless temperature and also its uncertainty: T* = (0.714 ±0.008). Therefore, the percentage uncertainty of this result is 1.1%. Also, the t-test calculated for this result showed Pð V=σ V Þ = 0.

Conclusion
The two-dimensional numerical model proposed for the cylindrical geometry in this article allowed determining simultaneously two thermal properties by optimization, their uncertainties and the correlation coefficient between them, using an experimental dataset of the cooling process of a cucumber. The results obtained for the average values of the parameters are compatible with results from the literature for the same product, and the statistical indicators chi-square and determination coefficient were χ 2 = 2.6981 ×10 −3 and R 2 = 0.9991, respectively. Due to the covariance matrix calculated during the optimization process, it was possible to simulate a new cooling curve for a similar cucumber, but with other dimensions, including a confidence band that allowed evaluating, by graphical inspection, the precision of this new simulation. Chi-square (dimensionless)