1.  Introduction

Nowadays, thermal sterilization of canned foods in retorts has been one of the most widely used methods of food preservation. Typically, the thermal sterilization method consists of heating food containers in pressurized retorts at defined temperatures for pre-established time intervals, followed by sudden cooling.

The conditions applied to the thermal processing of canned foods must be defined by specialists according to each set of equipment/container/food. This system defines the following process conditions: processing time, process temperature, minimal initial temperature of canned food, pressure profile, and deaeration program of the retort.

The determination of processing time is conducted with great care, according to the quality and safety requirements. This includes careful monitoring of processing time and temperature control to prevent under-processing or over-processing of canned food.

Although sterilization processes can be carefully defined, their day-to-day application in industrial production often deviates from ideal conditions. Operational problems, such as energy failure, temporary interruptions of steam supply, and failures of control system hardware, frequently occur and result in process deviations based on temperature variations of the heating medium. Other process deviations can occur when there are modifications in basic process parameters, such as the initial temperature of canned food, reduction of processing time, and lower heating temperature medium. When these deviations are observed, the batch is usually reprocessed or segregated for further analysis of the process trend records by specialists, who will determine the batch's fate. These procedures are costly and time-consuming, and when the decision is to reprocess the canned food, it implies a deterioration in the quality of the food product [1].

The sterilization of canned food partially depends on the heat penetration rate inside it. When food products are packaged into containers that will undergo heat treatment inside a retort, using steam as a heating agent, the condensation of the steam will determine the rate of heat transfer to the canned food/container system. The following factors affect the heat transfer rate: (i) heat transfer coefficients, (ii) physical properties of food products and containers, (iii) temperature difference between the retort and the canned food, and (iv) container dimensions [2].

Some information about the physical properties of canned foods is usually not available. Therefore, a simple method to estimate thermal diffusivity could be useful, as it influences the quality control of foods. [3] developed an experimental alternative method to determine the thermal diffusivity of wheat flour and cassava, considering them as models for analyzing a semi-infinite solid. The experiment involved an aluminum container filled with flour, an incandescent lamp as a heat source, and thermocouple sensors for acquiring temperature data. They also carried out a simulation in the Partial Differential Equation Toolbox of Matlab R2015a to compare the temperature distributions within the container.

A two-dimensional numerical model to determine thermal properties (their uncertainties and the covariance matrix) of cucumbers with cylindrical geometry, using an experimental dataset obtained during their cooling was proposed by [4]. The model is based on the Levenberg-Marquardt algorithm and can be used to determine parameters of a differential equation using experimental data.  The study provides insights into the behavior of food products during cooling and can be applied to simulate cooling curves for cucumbers under similar conditions, but with other dimensions, with no need for new experiments.

A simple, accurate, and robust method for determining the thermal diffusivity of solid foods was developed by [5]. The authors employed an inverse numerical method to estimate thermal diffusivity from transient temperature measurements at any position in a food sample, using a sequential parameter estimation technique. The developed algorithm was validated using both computer-generated temperatures and actual experimental temperature data. The estimated thermal diffusivity from the computer-generated data showed excellent agreement with the true value, and the estimated thermal diffusivity from the experimental measurements matched well with literature data.

According to [6], predicting the product temperature based on the temperature of the heating/cooling medium is the preferred approach in practice. However, it requires a mathematical model that accurately describes the temperature distribution in a transient regime for the food product. This approach offers advantages as it eliminates the need for installing test units to measure the product temperature, allowing for modeling and control of both continuous and batch processes. One significant disadvantage is the requirement of knowledge about the mode of heat transfer, as well as the relevant thermal and physical properties for each product and packaging type and geometry.

The studied product is a solid canned food, where heat conduction is the predominant heat transfer mechanism. The physical model is described by the heat equation.

                                                                                                           (1)

Where T is the food temperature and α is the thermal diffusion coefficient of the studied product which is a temperature function in this work.

Here, T represents the food temperature, and α is the thermal diffusion coefficient of the studied product, which is temperature-dependent in this work.

This study aims to provide reliable mathematical models for the sterilization process by identifying the thermal diffusivity of textured soy protein (TSP). This approach offers significant benefits to the field, ensuring safer consumption by reducing the risk of producing unsafe or overprocessed food. Additionally, it enables the production of high-quality food with lower energy consumption during sterilization and reduced nutrient degradation. The specific objectives of this study include determining the thermal diffusivity value of TSP during the sterilization process, using experimental data of retort temperature and the temperature in the geometric can center. The methodology employed involves applying a spatial discretization method based on finite volumes with high-order reconstruction in an unstructured mesh of triangles.

2.  Materials and methods

2.1. Differential equation resolution

An essential step in the parameter identification process is resolving the differential equation associated with the phenomenon (Equation 1). Several methods exist in the literature for solving partial differential equations. However, [7] state that most methods for solving PDEs typically offer a spatial precision order of two, which is widely used in academia and industries due to its intuitive nature, ease of implementation, and considerable precision for many flow problems. These methods also form the foundation for commercial fluid flow problem-solving software. However, when considering the parameter identification methodology, using numerical methods with a higher order of precision can reduce the error associated with the identified parameters.

 

The spatial discretization method employed in this work has several advantages. Firstly, it offers versatility since it can solve differential equations of any order, where the order n can be any natural number greater than or equal to three. Secondly, the method discretizes the domain using an unstructured mesh of triangles, making it more adaptable for equations defined in complex spatial geometries, such as solid revolution with azimuth symmetry. However, one disadvantage of this method is its implementation complexity. Using unstructured meshes requires the generation and storage of a complete data structure to identify the control volumes in the mesh [8].

Some authors have already applied the spatial discretization of PDEs using the finite volume method with high-order reconstruction in unstructured meshes. However, in the case of two-dimensional problems, only Cartesian coordinate systems were studied. In this work, an adaptation of this method is used for the cylindrical coordinate system, specifically in the radial and axial dimensions denoted by r and z, respectively.

The finite volume method with high-order reconstruction involves spatially describing the unknown functionT_i within the control volume  by employing a Taylor series expansion, while keeping the time t_j fixed (for j = 1, 2, ..., m). Here, m represents the number of points discretized in the temporal domain. This can be expressed by Equation 2 [9],

                         (2)

Where (r_i, z_i) represents the coordinates of the centroid of the i-th control volume (triangle),

  is the reconstructed solution at time t in the i-th control volume, and 

  denotes the partial derivatives of T evaluated at the centroid of the i-th control volume.

The approach used to obtain the reconstructed solution was adapted from [10, 11] to suit the cylindrical geometry.

An essential aspect of the finite volume method is the utilization of the integral formulation of the differential equation. Returning to Equation 1 and integrating both sides of the equation over the control volumeVi, we obtain:

                                                                                                                         (3)

Next, both sides of Equation 3 are divided by V_i , keeping in mind that , and the time derivative is interchanged with the volume integral, resulting in:

                                                                         (4)  

The term in parentheses in Equation 4 represents the mean value of temperature in control volume.  Since the reconstruction polynomial, and consequently the diffusion coefficient α, do not depend on the variable θ, it is possible to rewrite the second term using a contour integral through algebraic manipulations. This leads to Equation 5:

                                                                                                    (5)

wheredenotes the unit normal vector relative to the faces of each control volume.

Equation 5 represents a system of ordinary differential equations (ODEs) that can be expressed as:

                                                                                                                                                          (6)

  where .

The solution of the ODE system represented by Equation 6 requires a time integration methodology. Considering ease of implementation and lower computational cost per time step, an explicit third order Runge-Kutta method with four steps [12] was chosen. This method has the advantage of a larger stability region for the time step. The [12] method applied to the ODE system (Equation 6) is described by Equation 7:

                                                                             (7)

Since this is an explicit methodology, the method proposed by [12] (Equation 7) is numerically stable if the following condition is satisfied:

                                                                                                                                                     (8)

Where A_Li represents the area of the triangle corresponding to control volume i in the longitudinal section of the cylindrical domain.

2.2. Parameters identification methodology

According to [13], the thermal conductivity (κ) of several materials is better represented by a linear function than a constant value over a wide temperature range. For textured soy protein (TSP), it is reasonable to assume that the density and specific heat are constant [14]. Consequently, the thermal diffusivity value is also better represented by a linear function, as shown in Equation 9:

                                                                                                      (9) 

Where are the parameters to be identified.

Parameter identification problems are typically formulated as optimization problems. In the case of identifying the thermal diffusivity of TSP, the optimal values ofare those that make the model's temperature distribution closer to the experimentally obtained temperature distribution. Since only temperature values were collected at the geometric center of the can experimentally, the optimization problem involves minimizing a functional described by Equation 10:

                                                                                                              (10) 

where m is the number of points used, is the temperature value at the center of the can by simulation andis the value experimentally, that is, the reference value.

where m is the number of points used,represents the simulated temperature value at the can center, and  is the experimental value, i.e., the reference value.

 

The optimization problem was solved using the differential evolution method [15] with a population size of 10 individuals, a differential weight of F=0.8 and a crossover probability of CR=0.9 The DE/rand-to-best/1/Bin variant with parameter λ=0.5 [16] was employed.

 

As a stopping criterion, a stagnation condition was adopted, where the difference between two consecutive iterations of the objective function value (Equation 10) is less than10^-8 for ten consecutive iterations. To ensure more reliable results, the computational tests were performed in triplicate, and the best results are presented in this study.

All the numerical methods described above were implemented using the C++ programming language. The GNU GCC compiler version 13 was used, along with GNUplot software version 5.4 for generating graphs, and Gmsh software version 4.11 for generating the triangular mesh used.

 

2.3. Obtaining reference values

To obtain experimental data, [17] monitored the temperatures of TSP in a 300x407 can with a height of 0.12 m and a diameter of 0.072 m. They used a Data Trace temperature sensor placed at the geometric center of the can, which is the coldest point of the canned food during the sterilization process. This location is particularly susceptible to microbiological growth due to the lower thermal load at that point.

Due to the limitation imposed by the Spiteri and Ruuth method [12] for time integration (Equation 7), working with discrete data is not a good option. Therefore, the experimentally sampled temperature was adjusted using a first order plus dead time transfer function (FOPDT) [18], resulting in Equation 11. This adjusted temperature will be used as the reference for the minimization described in Equation 10.

                                              (11)

Furthermore, to optimize computational cost, only a portion of the data from the food heating step was used. These data were considered up to 2,000 seconds. Figure 1 shows a comparison between the experimental results for the heating step and their adjustments using the first order plus dead time transfer function.

 

Figure 1. Experimental results and FOPDT adjustments in the heating step of the TSP sterilization process.

 

1.  Results and discussion

In this section, we present the results obtained from the computational tests conducted to identify the parameters associated with the thermal diffusivity of TSP. A quarter of the cylindrical geometry was considered as the spatial domain due to the problem's symmetries. The initial condition used was the experimental temperature of 42.22 °C. The autoclave temperature, defined by Equation 12, was applied as a boundary condition for the region in direct contact with the steam.

                                                                 (12)      

At the axial axis of the can (r = 0), the temperature derivative was assumed to be zero since the temperature is minimal along that axis. For convenience, the upper right quarter of the can was used, and the z-axis scale was adjusted to z ∈ [0, 0.06] m.

All computational tests were performed using the developed C++ program. A mesh with 104 control volumes was utilized. [17] determined the thermal diffusion coefficient of TSP to be. This value served as a reference for defining the search space of the differential evolution method. The search space was divided for the analysis, considering one for the increasing values of the function α(T) and another for decreasing values.

Table 1 presents the parameters obtained for 3rd and 4th order simulations using the differential evolution method, along with the corresponding value of the objective function (Equation 10). The table also shows the value of the objective function for each order of spatial discretization.

Table 1. Parameters obtained from the optimization of Equation 10 compared to a constant value of α  

As expected, higher order simulations demonstrate lower errors, although the improvement is not significant enough to justify the use of orders greater than 4 due to the substantial increase in computational cost for a marginal precision gain.

The 2nd order simulation, which is the most common order, was not performed. The chosen method, using Gaussian quadrature, is unable to solve problems with constant initial conditions, such as this one. The derivatives of the reconstruction polynomial become constant values, resulting in a null residue R (T, t) in Equation 6. This does not occur for orders higher than 2.

Furthermore, it is observed that using a linear temperature function for the thermal diffusivity improves the fitting of the experimental data by the phenomenological model. The value of the objective function (Equation 10) is more than eight times smaller compared to the value obtained for the same order with constant thermal diffusivity.

Figure 2 displays the temperature values at the center of the can using the identified parameters. The results are compared with the simulation using constant diffusivity, the adjustment by FOPDT used as a reference, and the experimental data. The figure demonstrates the quality of the thermal diffusivity adjustment achieved by employing a linear temperature function.

 

Figure 2. Comparison between simulated temperatures, adjusted by FOPDT and experimental temperature for 3^rd and 4^th order.

 

In Figure 3, the transient temperature distribution of TSP is shown at different moments in time, using 4th order spatial precision with the parameters identified in Table 1. The 4th order precision was chosen for its slightly better results and a computational cost that is not significantly higher than the 3rd order. The transient distribution depicts the temperature evolution within the can over time, highlighting that the center requires the longest thermal processing time to reach the ideal temperature for ensuring food quality.

 

 

Figure 3. Simulation using parameters identified with 4th order in the Table 1 with reconstruction of also 4th order.

 

The analysis of the results obtained in this study provides a deeper understanding of the thermal diffusivity of textured soy protein (TSP) and its implications for the sterilization process of canned foods.

The identification of parameters related to thermal diffusivity enables more precise simulations and the ability to explore different scenarios and geometries by varying retort temperatures and heating times. This offers a valuable tool for optimizing sterilization processes, reducing costs and risks associated with physical testing, and ensuring high-quality and safe food for human consumption.

Furthermore, the improvement in fitting experimental data by using a linear function for thermal diffusivity highlights the importance of considering the temperature-dependent variation in diffusivity.

These findings can also provide insights for future research, such as incorporating other physical parameters and further experimental validation of the proposed models. In summary, the results of this study contribute significantly to the advancement of knowledge in the field of canned food sterilization and pave the way for new research and innovation opportunities.

 

4.  Conclusions

In this study, a numerical methodology based on the finite volume method with high-order reconstruction was developed to identify the thermal diffusivity parameters associated with textured soy protein (TSP). The methodology was implemented using a quarter of the cylindrical geometry due to the problem's symmetries, and a C++ program was developed for the computational tests.

The results obtained demonstrated the effectiveness of the proposed methodology in accurately estimating the thermal diffusivity of TSP. The optimization problem, formulated as a parameter identification problem, was successfully solved using the differential evolution method. The obtained parameters provided a temperature distribution that closely matched the experimental data, indicating the reliability of the identified thermal diffusivity values.

The numerical simulations revealed that higher-order precision led to improved accuracy, although the computational cost increased significantly. Therefore, the use of orders higher than 4 was deemed unnecessary for this problem. It was also observed that a linear temperature function for the thermal diffusivity resulted in a significantly better fit to the experimental data compared to a constant diffusivity.

The transient temperature distribution depicted in Figure 3 illustrated the evolution of temperature in the can over time, highlighting the extended processing time required for the center of the can to reach the desired temperature for food safety.

Overall, the developed methodology proved to be a valuable tool for the estimation of thermal diffusivity parameters in TSP. The accurate determination of these parameters is crucial for optimizing the sterilization process and ensuring the quality and safety of canned food products. Another applicability for the proposed model is the possibility of carrying out tests to change operating conditions, such as, for example, the temperature of the autoclave, without the need to carry out the physical experiment, which would save time and financial resources that would be spent to perform experimental tests. Future research can focus on applying the proposed methodology to other food materials and exploring additional optimization techniques to further enhance parameter identification accuracy.

Authors’ contributions

Conceptualization; investigation; software; validation; writing-original draft, R.Y.M.B.; Conceptualization; methodology; supervision; validation; writing-review and editing, A.A.S.; Conceptualization; investigation; methodology; software; writing-review and editing, R.G.

 

Acknowledgements

The authors gratefully acknowledge the financial support provided by the CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais).

 

Funding

Fundaçao de Amparo à Pesquisa do Estado de ˜Minas Gerais, Grant/Award Number: 02100-12/2013.

 

Availability of data and materials

All data will be made available on request according to the journal policy.

 

Conflicts of interest

Authors have declared that no competing interests exist.

 

 References

1.         da Silva, W.P.; Souto, L.M.d.; Ferreira, J. P.d.L.; Gomes, J.P.; Lima, A.G.B.d.; Queiroz, A.J.d.M.; Figueirêdo, R. M.F.d.; Santos, D.d.C.; Santana, M.d.F.S.d.; Santos, F.S.d. Use of inverse method to determine thermophysical properties of minimally processed carrots during chilling under natural convection. Foods. 2023, 12, 2084. https://doi.org/10.3390/ foods12102084.

2.         Teixeira, A.; Balaban, M.; Germer, S.; Sadahira, M.; Teixeira-Neto, R.; Vitali, A. Heat transfer model performance in simulation of process deviations. J. Food Sci. 1999, 64, 488–493. https://doi.org/10.1111/j.1365-2621.1999.tb15068.x.

3.         Mari, J.; Mari, M.; Ferreira, M.; Conceição, W.; Andrade, C. A simple method to estimate the thermal diffusivity of foods. J. Food Proc. Eng. 2018, 41, e12821. https://doi.org/10.1111/jfpe.12821.

4.         da Silva, W.P.; Silva e Silva, C.M.D.P.; de Souto, L.M.; Gomes, J.P.; Queiroz, A.J.M.; de Figueiredo, R.M.F. Thermal properties determination of a cylindrical product during its cooling: two-dimensional numerical model and uncertainty. Int. J. Food Prop. 2019, 22:1, 343-354. https://doi.org/10.1080/10942912.2019.1578791.

5.         Mohamed, I.O. Development of a simple and robust inverse method for determination of thermal diffusivity of solid foods. J. Food Eng. 2010, 101, 1– 7. https://doi.org/10.1016/j.jfoodeng.2010.05.002.

6.         Larousse, J.; Brown, B. Heat penetration in canned foods-theoretical considerations. Food Cann. Technol. 1997, 383–423.

7.         Vincent, P.E.; Jameson, A. Facilitating the adoption of unstructured high-order methods amongst a wider community of fluid dynamicists. Math.l Model. Nat. Phenom. 2011, 6, 97–140. https://doi.org/10.1051/mmnp/20116305.

8.         van Altena, M. High-order finite-volume discretisations for solving a modified advection-diffusion problem on unstructured triangular meshes. Master’s Thesis, University of British Columbia, Vancouver, Canada, 1999. https://doi.org/10.14288/1.0080945.

9.         Santana, A.A. Identificação de parâmetros em problemas de advecção-difusão combinando a técnica do operador adjunto e métodos de volumes finitos de alta ordem. Ph.D. Thesis. Universidade de São Paulo, São Paulo, Brazil, 2007. https://doi.org/10.11606/T.45.2007.tde-08012008-151101.

10.     Michalak, C., 2009. Efficient high-order accurate unstructured finite-volume algorithms for viscous and inviscid compressible flows. Ph.D. Thesis, University of British Columbia, Vancouver, Canada, 2009. https://doi.org/10.14288/1.0067108.

11.     Silva, J.E.C. p-Multigrid explícito para um método de volumes finitos de alta-ordem não estruturado. Ph.D. Thesis, Universidade de São Paulo, São Paulo, Brazil, 2016. https://doi.org/10.11606/T.45.2016.tde-24082016-174037.

12.     Spiteri, R.J.; Ruuth, S.J. A new class of optimal high-order strong-stability-preserving time discretization methods. SIAM J. Numer. Anals. 2002, 40, 469–491. https://doi.org/10.1137/S0036142901389025.

13.     Çengel, Y.A. Heat and Mass Transfer: A Practical Approach, 3rd ed., McGraw-Hill: New York, USA, 2006.

14.     Yano, T.; Kong, J.Y.; Miyawaki, O.; Nakamura, K. The “intrinsic” thermal conductivity of wet soy protein and its use in predicting the effective thermal conductivity of soybean curd. J. Food Sci. 1981, 46, 1357–1361. https://doi.org/10.1111/j.1365-2621.1981.tb04174.x.

15.     Storn, R.; Price, K. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization. 1997, 11, 341–359. https://doi.org/10.1023/A:1008202821328.

16.     Yang, X.S. Nature-inspired optimization algorithms, 1^st  Edn., Elsevier: London, UK, 2014.

17.     Gedraite, R.; Souza, D.L.; Kunigk, L. Simplified mathematical modeling of the heat transfer process in foods using the differential evolution algorithm, in: Lobato, F.S.; Steffen Jr, V.; Silva Neto, A.J. (Eds.), Computational Intelligence Techniques with Applications in Inverse Engineering Problems, 1st ed., OMNIPAX: Curitiba, Brazil, 2014, 11 – 22. https://doi.org/10.7436/2014.tica.0.

18.     Seborg, D.E.; Mellichamp, D.A.; Edgar, T.F.; Doyle III, F.J. Process dynamics and control, 3^rd Edn., John Wiley & Sons: Danvers, USA, 2010.