Use of Inverse Method to Determine Thermophysical Properties of Minimally Processed Carrots during Chilling under Natural Convection

The aim of this study was to determine the thermophysical properties and process parameters of cylindrical carrot pieces during their chilling. For this, the temperature of the central point of the product, initially at 19.9 °C, was recorded during chilling under natural convection, with the refrigerator air temperature maintained at 3.5 °C. A solver was created for the two-dimensional analytical solution of the heat conduction equation in cylindrical coordinates. This solver and the experimental data set were coupled to the LS Optimizer (V. 7.2) optimization software to simultaneously determine not only the values of thermal diffusivity (α) and heat transfer coefficient (hH), but also the uncertainties of these values. These values were consistent with those reported in the literature for carrots; in this study, the precision of these values and the confidence level of the results (95.4%) were also presented. Furthermore, the Biot numbers were greater than 0.1 and less than 40, indicating that the mathematical model presented in this study can be used to simultaneously estimate α and hH. A simulation of the chilling kinetics using the values obtained for α and hH showed good agreement with the experimental results, with a root mean square error RMSE = 9.651 × 10−3 and a chi-square χ2 = 4.378 × 10−3.


Introduction
The world production of carrot (Daucus carota L.) was about 42 million tons in 2021, concentrated mainly in China (18.2 million tons; 43.3% of world production) and, to a lesser extent, in Russia (1.3 million tons; 3.1% of world production) [1]. In a review study [2], the authors found that the sensory and chemical quality attributes of carrots can be influenced by several factors, such as genetic and climatic factors, cultivation method, storage temperature, and atmosphere, as well as thermal processing. Regarding the chemical composition of carrots, health-promoting bioactive compounds, such as vitamin C, phenolics (neochlorogenic and chlorogenic acid, 5-p-coumaroylquinic acid, or 3,5-dicaffeoylquinic acid), polyacetylenes (falcarinol, falcarindiol, and falcarindiol acetate), and, mainly, carotenoids (α-and β-carotene) [2][3][4][5][6][7][8], are worth mentioning. Conventionally, carrots are marketed whole, but the food industry has been dedicating a great deal of effort to developing products that combine safety, nutritional quality, and practicality in preparation. In this last aspect, minimally processed carrots are an attractive product option.
Minimal processing of fruits and/or vegetables consists of different steps, which include, but are not limited to, sanitization, removal of inedible parts, cutting or milling, packaging, and storage [9]. Despite the advantages related mainly to their easy preparation, minimally processed products (MPPs) are more perishable than their intact fresh versions due to the damage to their cellular tissue, which ends up being a window for microbial contamination and induces various deteriorating reactions [10][11][12], particularly under warmer storage conditions [13,14]. Thus, in order to extend the lifespan of MPPs, a fundamental part of the process is chilling the product and storing it under refrigeration [15][16][17][18][19]. In the chilling stage, products are typically at room temperature initially and are exposed to lower temperatures (from 2 to 10 • C) [20][21][22] until the thermal equilibrium condition is established.
Many phenomena involving food processing can be described through the diffusion equation. These phenomena include drying [23], osmotic dehydration [24,25], heating [26], and chilling [27]. To describe the chilling kinetics of a product, it is necessary to know its thermophysical properties, such as thermal diffusivity and the parameter heat transfer coefficient, among others. Thus, it is possible to estimate the minimum time for the least favorable point of the product to reach equilibrium temperature and, with this, determine the costs of the chilling process. These properties and parameters can be determined in two ways: directly and indirectly [28,29]. In the first one, direct, it is common to use specific equipment and relatively complex handling [30][31][32][33]. Conversely, the indirect way, based on the inverse method, has the advantage of being simpler and not requiring the use of any complex and/or expensive equipment. Briefly, the indirect way consists of solving, analytically or numerically, the partial differential equation of heat conduction and fitting the solution obtained, using some optimization process, to a set of experimental data of the temperature measured at a specific point of the product as a function of process time [34][35][36][37][38][39].
In a case of pasteurization, the difference between the temperature of the heated product and its initial temperature is relatively large (∆T > 40 • C) [40][41][42], and thermophysical properties should generally be considered variable [43]. In the case of chilling, the difference between the initial temperature and the temperature of the cooled product is relatively small (∆T < 20 • C) [27,[44][45][46]. Thus, in cases of chilling, in general, thermophysical parameters can be considered constants in the description of the transient process. Thus, the analytical solution of the heat conduction equation can be used to describe the chilling kinetics of agricultural products [45,47,48]. In these cases, the properties density (ρ), specific heat (C p ) and thermal conductivity (k) can be combined into only one, called thermal diffusivity, α = k/(ρC p ) [49,50], to describe the chilling kinetics. If the boundary condition of the heat conduction equation for the process is of the first kind, this single property α is sufficient to describe chilling. On the other hand, if the boundary condition appropriate to the problem is of the third kind, it will also be necessary to know the heat transfer coefficient [51].
It is interesting to highlight that the literature on the determination of the thermophysical properties of carrots, processed or not, is limited. However, one example is the study conducted by Chang and Toledo [52], who used the inverse method to estimate the thermal diffusivity and heat transfer coefficient during the heating of carrot cubes. In another example, Murakami [32] indirectly estimated the thermal diffusivity of cylindrical pieces of carrots subjected to different processing techniques. However, no studies were found reporting the determination of the thermophysical properties of minimally processed carrots (MPCs) during their chilling within a given temperature range. However, accu- rately determining thermophysical properties is important for designing and optimizing refrigerated storage systems for MPCs. In this context, this study aimed to determine the thermophysical properties (and their uncertainties) of cylindrical pieces of carrot during their chilling using the inverse method. Such determination involved: 1-a set of experimental data of temperature at the central point of the product during the transient process; 2-a solver developed for the analytical solution of the heat conduction equation; and 3-an optimization software called LS Optimizer for Differential Equations and Functions (http://www.labfit.net/LS.htm, accessed on 1 March 2023).

Materials
Carrots were acquired at the local market in the city of Campina Grande, Paraíba, Brazil. Only specimens with similar size and without physical damage or other injuries of any nature were washed in running water to remove surface dirt. Then, they were sanitized by immersion in chlorinated solution (200 ppm) for 15 min, rinsed in running water, and then dried with absorbent paper to remove surface moisture. Finally, the carrots were peeled and cut into cylindrical pieces (diameter D = 44.00 mm and length L = 40.00 mm) with the aid of a steel mold and using a stainless steel knife ( Figure 1). Diameter and length were measured using an INSIZE brand caliper, with a measuring range from 0 to 150 mm and precision of 0.01 mm.

Materials
Carrots were acquired at the local market in the city of Campina Grande, Paraíba, Brazil. Only specimens with similar size and without physical damage or other injuries of any nature were washed in running water to remove surface dirt. Then, they were sanitized by immersion in chlorinated solution (200 ppm) for 15 min, rinsed in running water, and then dried with absorbent paper to remove surface moisture. Finally, the carrots were peeled and cut into cylindrical pieces (diameter D = 44.00 mm and length L = 40.00 mm) with the aid of a steel mold and using a stainless steel knife ( Figure 1). Diameter and length were measured using an INSIZE brand caliper, with a measuring range from 0 to 150 mm and precision of 0.01 mm.

Chilling Experiment
The carrot pieces, initially at 19.9 °C, were subjected to chilling in a refrigerator (Brastemp ® ) at a temperature of 3.5 °C with natural air convection. Before placing a piece of carrot in the refrigerator chamber, the temperature of the device was adjusted and maintained for 24 h. It is important to mention that, to minimize possible temperature fluctuations, the interior of the refrigerator was previously lined with polystyrene and filled with previously cooled water bottles, already at the temperature of the experiment, T = 3.5 °C (Figure 2).

Chilling Experiment
The carrot pieces, initially at 19.9 • C, were subjected to chilling in a refrigerator (Brastemp ® ) at a temperature of 3.5 • C with natural air convection. Before placing a piece of carrot in the refrigerator chamber, the temperature of the device was adjusted and maintained for 24 h. It is important to mention that, to minimize possible temperature fluctuations, the interior of the refrigerator was previously lined with polystyrene and filled with previously cooled water bottles, already at the temperature of the experiment, T = 3.5 • C ( Figure 2).  Figure 2 presents an experimental setup that is useful on a laboratory scale for the determination of thermophysical properties but which obviously would not be recommended for industrial or commercial use for chilling purposes. However, the results obtained herein can be recommended, as will be observed below.
Automatic acquisition of temperature data at the central point of the cylindrical piece during chilling was performed using a system composed of a TH-095 Instrutherm digital thermometer, two K-type thermocouples, a USB-01 cable with an RS-232 serial port, and a Core i5 notebook. Thermocouple probes were inserted into the geometric centers of the cylindrical pieces of carrot, defined as r = 0 and y = 0, while a probe of another thermocouple was kept close to the carrot piece but in contact with the ambient air to record the temperature inside the refrigerator. Finally, the thermocouples were connected to the digital thermometer (TH-095), which sent (via USB cable) the temperature data to the notebook, which recorded these data in real time ( Figure 3).   Figure 2 presents an experimental setup that is useful on a laboratory scale for the determination of thermophysical properties but which obviously would not be recommended for industrial or commercial use for chilling purposes. However, the results obtained herein can be recommended, as will be observed below.
Automatic acquisition of temperature data at the central point of the cylindrical piece during chilling was performed using a system composed of a TH-095 Instrutherm digital thermometer, two K-type thermocouples, a USB-01 cable with an RS-232 serial port, and a Core i5 notebook. Thermocouple probes were inserted into the geometric centers of the cylindrical pieces of carrot, defined as r = 0 and y = 0, while a probe of another thermocouple was kept close to the carrot piece but in contact with the ambient air to record the temperature inside the refrigerator. Finally, the thermocouples were connected to the digital thermometer (TH-095), which sent (via USB cable) the temperature data to the notebook, which recorded these data in real time ( Figure 3).

Governing Equation
The partial differential equation describing heat conduction in a finite cylinder with radial and axial symmetry is given by and k is thermal conductivity (W m −1 K −1 ). In Equation (1), r is the radial position defined in relation to the central axis of the finite cylinder and, together with the axial position y (direction of the central axis), defines the location of a point inside the cylinder. Equation (1) can also be written as where α is thermal diffusivity (m 2 s −1 ) and is given by k/(ρC p ). thermometer, two K-type thermocouples, a USB-01 cable with an RS-232 serial port, and a Core i5 notebook. Thermocouple probes were inserted into the geometric centers of the cylindrical pieces of carrot, defined as r = 0 and y = 0, while a probe of another thermocouple was kept close to the carrot piece but in contact with the ambient air to record the temperature inside the refrigerator. Finally, the thermocouples were connected to the digital thermometer (TH-095), which sent (via USB cable) the temperature data to the notebook, which recorded these data in real time ( Figure 3).  To solve Equation (2) analytically, the following assumptions were made: (1) Heat generation within the product during chilling can be considered negligible; (2) initial temperature is known and uniform; (3) the geometry of the product is that of a homogeneous and isotropic finite cylinder; (4) the origin of the coordinate system is located in the geometric center of the product; (5) for the chilling temperature range, the thermophysical properties of the product can be considered constant; (6) the appropriate boundary condition for the heat conduction equation is of the third kind, and the heat transfer coefficient is the same for all product/air interfaces; and (7) the chilling temperature T eq can be considered constant throughout the process, with cold air under natural convection.

Initial Condition and Radial and Axial Symmetry
For t = 0, where T 0 is the initial temperature (K), R is the cylinder radius (m), and L is the cylinder length (m). On the other hand, in the solution of the heat conduction equation, the assumptions of radial and axial symmetry were considered.

Boundary Condition
The boundary condition of the third kind presupposes resistance to heat flow at the boundary of the product. Thus, the equality between the internal heat flow (diffusive) on the surface of the product and the external flow (convective) in the outer vicinity of this surface is given by Equation (4) for the radial case (r = R): h H is the heat transfer coefficient (W m −2 K −1 ). For axial flows (y = ±L/2), the boundary condition of the third kind is given by Equation (5): where L is cylinder length (m) and T(r, y, t)| y=±L/2 is temperature (K) in the circular areas located in positions y = ±L/2.

Analytical Solution of the Heat Conduction Equation
Considering the simplifying assumptions presented above, the solution of Equation (2) can be obtained by the separation of variables [53,54].
where T(r,y,t), T 0 , and T eq are: (1) the temperature in a position (r,y) of the cylinder at an instant t; (2) the initial temperature of the product; and (3) the equilibrium temperature after chilling. In Equation (6), the parameter α is thermal diffusivity, and the coefficients A n and A m are calculated as follows: and In Equations (7) and (8), Bi 1 and Bi 2 are the Biot numbers for heat transfer in an infinite cylinder and on an infinite wall, respectively, and are calculated according to Equations (9) and (10), in this order: and In Equations (9) and (10), h H is the heat transfer coefficient, which was considered to have the same value for all cylinder surfaces. In Equations (6)-(8), µ n and µ m are the roots of the characteristic equations for an infinite cylinder and an infinite wall, respectively, and were calculated using the following equations: and In Equation (11), J 0 and J 1 are the first-kind Bessel functions of orders 0 and 1, respectively. On the other hand, Equation (6) can be defined in terms of a dimensionless temperature as follows:  (6) is used to determine the temperature T(r,y,t), whereas Equation (13) can also be used to determine the dimensionless temperature T*(r,y,t), for previously defined values of h and α. For this, the Biot numbers were calculated by Equations (9) and (10). The first 200 roots of Equations (11) and (12) were calculated by the bisection method, with a convergence tolerance of 10 −12 . This allowed determining the first 200 coefficients, A n and A m , given by Equations (7) and (8), respectively. Thus, for the analytical solution of the partial differential equation of heat conduction, a solver was created in FORTRAN. The developed software was created in Compaq Visual Fortran Professional Studio, V. 6.6.0, using a programming language option called QuickWin Application.

Inverse Problem: Simultaneous Determination of Optimal Values of α and h H
The LS Optimizer software, developed by the first author of this article, was used to determine the values of α and h H . For this, the solver created according to the previous sections and the experimental data set of temperature at the central point of the product were coupled to the LS environment. It is worth pointing out that LS Optimizer is ready-to-use software that uses non-linear regression (least squares) to determine parameters of an ordinary or partial differential equation (and functions) through known experimental data, using the Levenberg-Marquardt algorithm [55,56]. LS runs the solver to obtain the information needed for the determination of parameters. Thus, starting from initial values α 0 and h 0 , LS provides the optimal values for α and h H , their uncertainties, and the covariance matrix. In this study, the initial values of α and h were established as α 0 = 1 × 10 −7 m 2 s −1 and h H0 = 4 W m −2 K −1 . It is interesting to mention that the solver to be run by LS Optimizer was developed to read an "exp.txt" file, which contains the experimental data in three columns: the first with data for the independent variable (t), the second with values for the dependent variable ( T exp i at the central point, and the third column with the uncertainty values of the dependent variable (σ T i ), if known. However, if this last piece of information is not known, the user just needs to write the value equal to 1, and the LS already identifies the absence of these uncertainties. It is worth mentioning that the developed solver must read a file called "parameters.txt", generated by LS Optimizer, with the information given below. The first piece of information is an integer variable called "information", which can be, in this case, 0, 1, or 2. The specific value defines the name of the file to be generated by the solver, which can be unsteady.txt (if the value is 0), unsteady_a1.txt (if the value is 1), or unsteady_a2.txt (if the value is 2). Thus, after the information provided by the LS is identified, the solver stores the simulation results in the appropriate file for subsequent calculations of the sensitivities of the parameters to be determined. The second piece of information is an integer variable called "N_Param" that indicates the number of parameters involved in the solution of the differential equation (in this case, N_Param = 2). Finally, the third piece of information indicates the values established by the LS Optimizer for the parameters in each iteration (in this case, only two parameters: a 1 = α and a 2 = h). Thus, the solver calculates the dependent variable for the same values as the independent variable specified in the experimental dataset ("exp.txt") and then writes the simulations obtained in the files mentioned above (unsteady.txt, un-steady_a1.txt, or unsteady_a2.txt). These files are then used by LS Optimizer to calculate the chi-square in each iteration χ 2 = ∑ T i exp − T i sim 2 /σ T i 2 and also the corrections (∆a 1 and ∆a 2 ) of the fitting parameters (a 1 and a 2 , respectively). With the corrections to the parameters determined, the new values of the fitting parameters a 1 (a 1 = a 10 + ∆a 1 ) and a 2 (a 2 = a 20 + ∆a 2 ) are calculated. These corrected values are then considered new initial values (a 10 and a 20 ), and the iterative process occurs until a convergence criterion (|∆a 1 /a 1 | < 1 × 10 −5 and |∆a 2 /a 2 | < 1 × 10 −5 ) is reached. As additional information, the LS Optimizer software is freely available for download at the link provided at the end of the introduction to this article. On the other hand, in each iteration, the corrections ∆a 1 and ∆a 2 are calculated by minimizing the objective function (in the present case, chi-square, χ 2 ) [44], resulting in the following system of equations given in matrix form: where and The sums range from i equal to 1 to N, which is the number of experimental points. Additionally, (17) where "i" is the ith experimental point. In the last equations, T and Silva et al. [43] can be consulted. As additional information, when the relative tolerance on the two parameters is reached in the iterative process, the inverse of Equation (15) is the covariance matrix.

Results
Several tests were carried out with different pieces of carrot, whose results agreed with each other. However, we have chosen only one dataset for the analyses: the one with the least noise during data acquisition. Thus, the acquisition of temperature data at the central point of the cylindrical piece of carrot under natural convection of cold air resulted in the graph shown in Figure 4, which contains 49 experimental points.
At this point, it is important to emphasize that the values of the uncertainties of the properties and parameters, objects of this study, are not determined by repetitions of the measurement process (replicates), but through the covariance matrix, [M] −1 . This matrix is the inverse of the last matrix [M] calculated by Equation (15) in the iterative process, that is, the one obtained in the convergence of the iterative process referring to the non-linear regression.
The temperature data throughout the transient regime were then written in the dimensionless form, as explained by Equation (13), T * (r, y, t)= (T (r, y, t) − T eq )/(T 0 − T eq . Thus, the results obtained according to the proposed methodology were presented and analyzed as shown below.

Results
Several tests were carried out with different pieces of carrot, whose results agreed with each other. However, we have chosen only one dataset for the analyses: the one with the least noise during data acquisition. Thus, the acquisition of temperature data at the central point of the cylindrical piece of carrot under natural convection of cold air resulted in the graph shown in Figure 4, which contains 49 experimental points.

Thermophysical Properties and Parameters
For the chilling experiment, the initial temperature of the cylindrical piece of carrot was 19.9 • C, and the chilling temperature of the refrigerator chamber was maintained at approximately 3.5 • C. The cylindrical piece of carrot had a mass of 61.0 g and a moisture content on a wet basis of X wb = 90.4%. Its radius was R = 0.022 m, and its length was L = 0.040 m. Thus, the density could be calculated directly from the equation ρ = m/V, resulting in ρ = 1003 kg m −3 . The specific heat of the product was estimated by the model of Fikiin [58], C p = 1382 + 2805 X wb /100, and the result was C p = 3918 J kg −1 K −1 . Table 1 shows α and h, determined by optimization using: (1) the experimental data set of dimensionless temperatures at the center of the product during its chilling; (2) the solver developed to predict these dimensionless temperatures for α and h stipulated by the LS Optimizer; and (3) the LS Optimizer optimization software. The values obtained for α and h were multiplied by ρC p to obtain the thermal conductivity, k, and the heat transfer coefficient, h H , respectively. Such results are also presented in Table 1. The result obtained through the proposed model for the thermal diffusivity value of carrot was α = (1.43 ± 0.18) × 10 −7 m 2 s −1 , with a confidence level of 95.4%. This result was consistent with the value obtained from the estimate by Riedel correlation (α = 0.88 × 10 −7 + 0.6 × 10 −7 × X wb /100) [59], α = 1.42 × 10 −7 m 2 s −1 , and also with the value that Murakami [32] reported for fresh carrots (X wb = 87.6%; 25 • C): α = (1.44 ± 0.01) × 10 −7 m 2 s −1 . Chang and Toledo [52] determined the thermal diffusivity, α = (1.86 ± 0.04) × 10 −7 m 2 s −1 , of carrot pieces with a moisture content of X wb = 89.2%, at 138 • C. However, this estimate, at such a high temperature, could not really have good agreement with the results obtained in the present study, whose temperature range was 19.9 to 3.5 • C. On the other hand, the thermal conductivity value obtained using the sweat correlation [60] (k = 0.148 + 0.493 X wb /100) was k = 0.594 W m −1 K −1 . This value is compatible with that obtained in the present study (k = αρC p ), which was k = (0.56 ± 0.07) W m −1 K −1 , with a confidence level of 95.4%.
It is interesting to point out that, although Murakami [32] determined the uncertainty of the α value by error propagation (considering that the properties k, ρ, and C p are determined independently of each other, that is, covariance equal to zero between each pair of values), the confidence level was not specified. In the case of the study conducted by Chang and Toledo [52], the uncertainty of the result obtained for α was determined using 10 repetitions, which implies a confidence level of less than 68%. On the other hand, in the present study, it was possible to determine the values of α and h (and their uncertainties) using 49 experimental points obtained for the dimensionless temperature of the center point of the cylindrical piece during the transient process. In addition, the uncertainties originally provided by the LS Optimizer program for α and h were multiplied by the factor 2, which was recommended by the optimization software to ensure the 95.4% confidence level (assuming a Gaussian distribution of errors). It should also be noted that an additional advantage of the proposed model is the possibility of determining the values of Biot numbers for the infinite cylinder (Bi 1 ) and infinite wall (Bi 2 ) (Table 1), referring to the chilling of carrot pieces with the geometry of a finite cylinder. The low values of Biot numbers (Bi 1 = 0.271 and Bi 2 = 0.246), presented in Table 1, indicate the importance of considering the external resistance to heat transfer between the product and the cold air around its vicinities. Thus, it can be inferred that the boundary condition of the third kind, used in this study, is actually the appropriate boundary condition to describe the chilling process of cylindrical pieces of carrot under the natural convection of cold air. Table 1. Thermophysical properties and parameters determined for a carrot piece (X wb = 90.4%) during chilling from 19.9 • C down to 3.5 • C. The intervals are given with a 95.4% confidence level.

Properties
Carrot Reference The covariance matrix provided by the LS program is given by Equation (18).
[M] −1 = 8.2578 × 10 −17 −1.8611 × 10 −16 −1.8611 × 10 −16 5.2435 × 10 −16 (18) The covariance matrix, given by Equation (18), indicates that the parameters α and h are dependent on each other, and this dependence results in a correlation equal to −0.8944. This means that when the value of one of the parameters tends to increase, the value of the other parameter tends to decrease, which is common in optimization procedures involving two parameters. In this case, in error propagation calculations involving α and h, the covariance between these parameters must necessarily be considered. Similar results for correlation were found in several studies involving the determination of two parameters that describe phenomena involving the processing of food products [23,45,61,62].

Product Chilling Kinetics
With the results obtained for α and h, it was possible to simulate the chilling kinetics for the central point of the carrot piece, as shown in Figure 5.
A visual inspection of Figure 5 shows that the fit of the simulated line to the experimental data can be considered good. This qualitative observation is also verified by the high value of the determination coefficient (R 2 = 0.9991), the low values of chi-square (χ 2 = 4.378 × 10 −3 ), root mean square error (RMSE = 9.651 × 10 −3 ), and the standard deviation of the fit (σ = 9.25 × 10 −3 ). On the other hand, the determination of the values for α and h H allows simulating the chilling kinetics not only in the center of the carrot piece, at coordinates (0,0), but also at any other point, such as at the boundary of the product, at the coordinate point (22,0). A graph that shows these two simulations is presented in Figure 6. This figure shows what is expected: in intermediate times of chilling kinetics, the dimensionless temperature of the center is greater than the dimensionless temperature on the surface of the product, and this difference was maximum (∆T* = 0.11) at the instant t = 633.6 s (that is, t = 10.56 min).

Product Chilling Kinetics
With the results obtained for α and h, it was possible to simulate the chilling kinetics for the central point of the carrot piece, as shown in Figure 5.  A visual inspection of Figure 5 shows that the fit of the simulated line to the experimental data can be considered good. This qualitative observation is also verified by the high value of the determination coefficient (R 2 = 0.9991), the low values of chi-square (χ 2 = 4.378 × 10 −3 ), root mean square error (RMSE = 9.651 × 10 −3 ), and the standard deviation of the fit (σ = 9.25 × 10 −3 ). On the other hand, the determination of the values for α and ℎ allows simulating the chilling kinetics not only in the center of the carrot piece, at coordinates (0,0), but also at any other point, such as at the boundary of the product, at the coordinate point (22,0). A graph that shows these two simulations is presented in Figure 6. This figure shows what is expected: in intermediate times of chilling kinetics, the dimensionless temperature of the center is greater than the dimensionless temperature on the surface of the product, and this difference was maximum (ΔT* = 0.11) at the instant t = 633.6 s (that is, t = 10.56 min). It is interesting to highlight that, although the results observed in Figure 6 are expected, the proposed model can predict in detail such results for this transient process. It is also interesting to observe in Figure 6 that the temperature on the surfaces of the cylinder decreases slowly, a boundary condition of the third kind for the description of carrot chilling under natural convection. Thus, low Biot numbers could have been expected, as indeed they were.
The solver developed for the proposed model also allowed predicting temperature distributions in a circle at a given axial position y within the finite cylinder representing the carrot piece, which is shown in Figure 7a. Such distributions can be seen in Figure 7b for circles in different positions between y = −20.00 mm and y = +20.00 mm at time t = 13 min. For y = 0, for example, it is possible to observe that the chilling front advances from the outer surface toward the center of the product. It is interesting to highlight that, although the results observed in Figure 6 are expected, the proposed model can predict in detail such results for this transient process. It is also interesting to observe in Figure 6 that the temperature on the surfaces of the cylinder decreases slowly, a boundary condition of the third kind for the description of carrot chilling under natural convection. Thus, low Biot numbers could have been expected, as indeed they were.
The solver developed for the proposed model also allowed predicting temperature distributions in a circle at a given axial position y within the finite cylinder representing the  Figure 7a. Such distributions can be seen in Figure 7b for circles in different positions between y = −20.00 mm and y = +20.00 mm at time t = 13 min. For y = 0, for example, it is possible to observe that the chilling front advances from the outer surface toward the center of the product. This type of temperature distribution has also been described in studies on the chilling of different agricultural products, such as fig [27], cucumber [63], banana [45], potato [64], and apple [46]. Although this is a fairly obvious phenomenon, the proposed model makes it possible to predict these dimensionless temperature distributions at any given instant, as can be seen, for example, for the circular area in y = 0, during the first 13.0 min (780 s) of chilling, through the following link: http://labfit.net/Cooling.gif (accessed on 1 March 2023). Once again, it is evident from Figure 7b that, indeed, the boundary condition appropriate for the description of carrot chilling under natural convection is the third kind. Additionally, the animation shown via the link above clearly indicates that there are temperature gradients inside the carrot piece, with the temperature at the boundary varying over time.
According to many authors, such as Brischetto [65], if Bi < 0.1, a lumped model can be used to describe heat transfer problems with these Biot numbers, and the error associated with using this model to describe such problems is small. In this case, only the heat transfer coefficient can be determined. On the other hand, if Bi > 40, which corresponds practically to an infinite value, the appropriate boundary condition is of the first kind (and only the property of thermal diffusivity can be determined).
Unlike natural convection, it is interesting to observe that, depending mainly on the velocity of the coolant fluid, the heat transfer can be better described by the diffusion equation with the first kind of boundary condition. In this case, the equations describing the chilling process are the same, Equation (6) or (13), but the coefficients given by Equations (7) and (8) Additionally, when Bi1 and Bi2 tend to be infinite (in practical terms, Bi > 40), Equations (11) and (12) become J 0 µ n = 0 and cot µ m = 0, respectively.

Conclusions
As Bi1 and Bi2 are greater than 0.1 and both values are less than 40, it can be concluded This type of temperature distribution has also been described in studies on the chilling of different agricultural products, such as fig [27], cucumber [63], banana [45], potato [64], and apple [46]. Although this is a fairly obvious phenomenon, the proposed model makes it possible to predict these dimensionless temperature distributions at any given instant, as can be seen, for example, for the circular area in y = 0, during the first 13.0 min (780 s) of chilling, through the following link: http://labfit.net/Cooling.gif (accessed on 1 March 2023). Once again, it is evident from Figure 7b that, indeed, the boundary condition appropriate for the description of carrot chilling under natural convection is the third kind. Additionally, the animation shown via the link above clearly indicates that there are temperature gradients inside the carrot piece, with the temperature at the boundary varying over time.
According to many authors, such as Brischetto [65], if Bi < 0.1, a lumped model can be used to describe heat transfer problems with these Biot numbers, and the error associated with using this model to describe such problems is small. In this case, only the heat transfer coefficient can be determined. On the other hand, if Bi > 40, which corresponds practically to an infinite value, the appropriate boundary condition is of the first kind (and only the property of thermal diffusivity can be determined).
Unlike natural convection, it is interesting to observe that, depending mainly on the velocity of the coolant fluid, the heat transfer can be better described by the diffusion equation with the first kind of boundary condition. In this case, the equations describing the chilling process are the same, Equation (6) or (13), but the coefficients given by Equations (7) and (8) are defined as [44]: and Additionally, when Bi 1 and Bi 2 tend to be infinite (in practical terms, Bi > 40), Equations (11) and (12) become J 0 (µ n )= 0 and cot µ m = 0, respectively.

Conclusions
As Bi 1 and Bi 2 are greater than 0.1 and both values are less than 40, it can be concluded that the mathematical modeling presented in this article is adequate to describe the chilling of the carrot pieces under natural convection. On the other hand, the thermophysical properties of cylindrical carrot pieces were determined using non-linear regression during a transient process. The result obtained for the thermal diffusivity value of the product was α = 1.43 × 10 −7 m 2 s −1 , using the proposed model: (1) Data set of temperature in the center of the product; (2) solver developed; (3) LS Optimizer software. This result is consistent with the values found in the literature for similar temperatures of the product. However, the advantage of the proposed model is that it provides the temperature range within which this value is valid (from 19.9 down to 3.5 • C) and, mainly, its uncertainty with a confidence level of 95.4%. From the results obtained for the heat transfer coefficient (h H = 6.92 W m −2 K −1 ) and also for the Biot numbers (Bi 1 = 0.272 and Bi 2 = 0.246), it can be concluded that the surface resistance to heat flow should be considered in the description of the chilling process of carrot pieces with natural convection of cold air. Thus, the proposed model was adequate to describe the chilling of carrot pieces with the geometry of a finite cylinder, and the boundary condition appropriate to the process is of the third kind.