Numerical simulation of mushrooms during freezing using the FEM and an enthalpy: Kirchhoff formulation

The shelf life of mushrooms is very limited since they are susceptible to physical and microbial attack; therefore they are usually blanched and immediately frozen for commercial purposes. The aim of this work was to develop a numerical model using the ﬁnite element technique to predict freezing times of mushrooms considering the actual shape of the product. The original heat transfer equation was reformulated using a combined enthalpy-Kirchhoff formulation, therefore an own computational program using Matlab 6.5 (MathWorks, Natick, Massachusetts) was developed, considering the difﬁculties encountered when simulating this non-linear problem in commercial softwares. Digital images were used to generate the irregular contour and the domain discretization


List of symbols CG
Global capacitance matrix Cp Specific heat J (kg °C) -1 Cp ap Apparent specific heat J (kg °C) - The shelf life of mushrooms, such as the commercial button mushroom Agaricus bisporus, is limited to a few days, mainly because they have no cuticle to protect them from physical or microbial attack and water loss.In the same way, they have a large content in nutrients and high respiration rate; factors that induce deterioration immediately after harvest [1].Finally their high tyrosinase and phenolic content makes them very susceptible to enzymatic browning [2].In view of their highly perishable nature, mushrooms must be processed to extend their commercial shelf life for off-season use [3].The conventional conservation process of mushrooms includes washing of the samples and immediate blanching in order to inactivate enzymes, induce a volume contraction, and make the product more pliable to facilitate filling operations [4].During this step shrinking of the tissue samples reaches up to 18%, as reported in [5].Then mushrooms can be submitted either whole or sliced to several processes in order to obtain longer durability, for example low (refrigeration or freezing) or high (sterilization, pasteurization, dehydration, etc.) temperature thermal processes.The process of freezing is generally regarded as being superior to canning or dehydratation when judged on the basis of retention of sensory attributes and nutritive properties [6], and it combines the favorable effect of low temperatures with the conversion of water into ice [7].Food engineers are interested in predicting cooling and freezing times in order to estimate the energy requirements for freezing and to design the necessary equipment for effective and rational processing.For this purpose, the evolution and distribution temperatures in the whole dominium of the food during freezing must be known.To this end, finite element method (FEM) is an established formulation for numerical predictions of the transient temperature in heat conduction problems [8] such as chilling and freezing.It has the great advantage that it can deal better with problems where the object has an irregular geometry [9,10].On the other hand, the most important issue regarding freezing tunnels is the knowledge of the coefficients of heat and mass transfer between product and air, which are necessary for the design of industrial equipment and the heat transfer operations [11].During the freezing process, which involves the phase change of water into ice in the food product, the thermophysical properties such as specific heat, thermal conductivity and density undergo abrupt changes due to the latent heat release.The system is then established as a highly non-linear mathematical problem.
Several techniques were applied to deal with the large latent heat release when using the FEM.One of the traditional methods is the use of the apparent specific heat, where the sensible heat is merged with the latent heat to produce a specific heat curve with a large peak around the freezing point, that can be considered a quasi delta-Dirac function with temperature (depending on the amount of water in the food product).The abrupt change in the apparent specific heat curve requires several iterations for each time step and usually destabilizes the numerical solution.In some commercial simulation softwares that implement the FEM, such as COMSOL Multiphysics, the use of the apparent specific heat is the only method available [10].Many authors have done approximations by ''softening'' the peak curve in order to obtain some convergence of the method, modifying the shape of the apparent specific heat curve while maintaining the total latent heat constant.However, this softening method is not recommended, because the actual temperature range around the freezing zone, is altered becoming wider than the actual temperature freezing range.
Sheen and Hayakawa [12] have simulated freezing of whole mushroom shapes using Tylose gel material with 77% water content.Numerical problems were also reported in this work when applying the finite difference numerical technique caused by an unrealistic heat balance around each surface node on the curvilinear boundary.
Sliced mushrooms which have an irregular 3D geometry are also submitted to freezing, however this type of shape were not considered in previous works.
The implementation of the enthalpy method, which can be obtained through the integration of the specific heat with temperature [10,13,14] and the Kirchhoff function, which is the integral of the thermal conductivity, allows the reformulation of the heat transfer differential equation into a transformed partial differential system with two mutually related dependent variables H (enthalpy) and E (Kirchhoff function) [15][16][17][18].Even though it generates great advantages to the resolution of the phase change problem, there are few works in literature which apply this method to simulate either freezing or thawing of other foodstuffs.
Combining both transformations helps to avoid inaccuracies and/or divergence of the numerical method, caused by the latent heat peak release and the jump of the thermal conductivity at the phase transition, with the great numerical advantage of minimizing the execution speed of the program, since the resulting finite element matrices are constant.
The goals of this work are: 1. To develop a finite element code that simulates: the freezing process of an irregular shaped food using a combined enthalpy and Kirchhoff transformation method.2. To apply the numerical model during freezing of whole and sliced mushrooms considering a twodimensional axial symmetric domain and a 3D geometry, respectively.3. To validate the numerical solutions comparing the temperature predictions with experimental data during the freezing process of whole and sliced mushrooms in a tunnel.4. To predict processing times for different operating industrial conditions; different external fluid temperature and surface heat transfer coefficients considering both shapes.
2 Materials and methods 2.1 Finite element simulation during the freezing process of a 3D geometry (mushroom slice) During the phase change transition in the freezing process the thermo-physical properties are strongly dependent on temperature.This constitutes a highly non-linear mathematical problem.A mushroom slice represents an irregular three dimensional geometry, therefore the heat conduction equation in Cartesian coordinates with phase change transition can be written as follows: Equation 1 is valid in the domain X, where T is the temperature, k is the thermal conductivity, Cp the apparent specific heat, and q the density [19].The initial (Eq.2) and boundary conditions (Eq. 3) are: where dX is the domain of the convective interface which corresponds to external surface in contact with the cooling air, n x , n y , and n z are the normal outward vector components, T ext is the external air temperature, T 0 is the initial food temperature and h is the surface heat transfer coefficient.
By performing the following change of variables: where H is defined as the volumetric specific enthalpy [20], E is the Kirchhoff function that represents the thermal conductivity integral [15,20], and T* is a reference temperature that corresponds to a zero value of H and E. By combining Eqs. 1, 4 and 5, and the initial and boundary conditions represented by Eqs. 2 and 3 the following equations were obtained: Applying the finite element technique the following system of ordinary differential equations must be solved: where: ðN T NÞdX e is the global capacitance matrix ðB T BÞdX e is the global conductance matrix H, E, and T are the nodal values of enthalpy, the Kirchhoff function, and temperature, respectively.N is the vector of dimensions [1 9 4] containing the shape functions (N j ) with j = 1-4 for the reference tetrahedron element, N t is the transpose vector (dimension 4 9 1), e is the total number of elements, e1 is the number of boundary elements, X e is the integration domain, dX s is the boundary integration domain.The matrix B (dimension 3 9 4) is defined as follows: It can be observed that the thermal properties of the foodstuffs have been incorporated into the new mathematical formulation (Eqs.6-8) as dependent variables; therefore these matrices need to be calculated only once, which reduces the computer time significantly.The extraordinary advantage of the Enthalpy method is that instead of solving the original partial differential equations using a temperature variable the problem was reformulated into another partial differential equation which uses an Enthalpy variable, and when implementing the variational formulation using the FEM the problem to be solved can be represented as a System of Ordinary differential Equations.Therefore the solution given as enthalpy values over the entire domain is a function H = H(time, x, y, z) If the surface heat transfer coefficient is considered constant during the process, FG and m also remain constant.
Equation 9 is also a system of ordinary differential equations with three unknown variables, H, E, and T that are interrelated through non-linear algebraic functions (H(T), E(T), H(E), T(H), T(E), E(H)).The temperature dependence of the thermal properties are incorporated into the enthalpy formulation by using this T(H) and E(H) functions.First the H(T) and E(H) were constructed according to Eqs. 4, 5 with the Cp(T) and k(T) integrals.The estimation of the temperature dependence thermal properties in mushrooms, Cp(T) and k(T), is explained in the following Section.The H(T) and E(H) functions are positive and monotonically increasing for each value of H there is only one value of T and E, i.e. they are bijective functions.Next using Matlab 6.5 interpolating functions the T(H) and E(H) functions were obtained.The functions were inserted into the ODE through the interp command in Matlab language.
Incorporating the functions E(H) and T(H) in Eq. 9 the system can be rewritten as: The system was solved using the standard Matlab 6.5 routines ODE (Ordinary Differential Equations) were the enthalpy values at each node are calculated for each time step.The temperature distribution was obtained using the function T(H), previously determined.As a result the temperature over the domain can be estimated using the interpolating function used in the finite element technique; this constitutes a great advantage over the finite difference where the temperature is given only at single points in the domain.

Finite element simulation during the freezing process of a 2D axial symmetric geometry (whole mushroom)
A whole mushroom can be approximated as a solid of revolution having a symmetry plane at r = 0. Therefore using cylindrical coordinates the equation that constitutes the heat transfer with phase change transition is: The boundary and initial conditions are: Using the Enthalpy and Kirchhoff transformation and afterwards the finite element variational formulation the same semi-discrete system as Eq. 9 is obtained: where : H, E, and T are the nodal values of enthalpy, the Kirchhoff function, and temperature, respectively.The matrix B (dimension 2 9 3) is defined as follows: This system was solved using the standard Matlab 6.5 routines ODE (Ordinary Differential Equations) which minimizes the computational efforts of the numerical algorithm.

Reconstruction of mushroom shape and mesh generation
Two-dimensional axial symmetric domain (whole mushrooms) and a 3D geometry (mushroom slices) were built from images of samples (Figs.1a and 2a).These images were digitally processed to obtain a binary image.In this regard, Image Processing Toolbox, MATLAB 6.5 (Math-Works, Natick, Massachusetts), was used according with the following steps [21]: -Conversion of original RGB images to grey-scale format, -Noise reduction through a 3 9 3 median filter to enhance image quality, -Segmentation through a threshold value which was obtained by analyzing the grey-scale image histogram.
A binary image was obtained where black colour (pixel value equal to 0) represented the background and white colour the sample (pixel value equal to 1).
The contour of the binary image was approximated with a B-Spline curve, which was used as a base for the construction of the simulation domain.To construct the axialsymmetric two-dimensional domain, the B-Spline curve was transformed into a solid object.The three-dimensional domain was obtained by extrusion of the two-dimensional irregular image.

Mesh details
To run the finite element model the two-dimensional axial-symmetric domain and the three-dimensional domain were imported into a mesh generator and discretized using triangles and tetrahedrons, respectively (Figs. 1b and 2b).
As shown in these figures a non-uniform grid system was used in the simulations.An unstructured mesh with 501 nodes and 908 triangular elements was developed for 2D model, while for the 3D model a mesh with 1,882 nodes and 7,693 tetrahedral elements was used.To achieve this meshing, a maximum element size of 1 mm and an element growth rate of 1.5 were specified for both cases.This will give the adequate number of elements.The use of finer mesh showed no significant effect on the accuracy of the solution.
The mesh information was preprocessed and used as input into the main finite element code.

Thermo-physical properties of the product
The numerical modelling of the freezing process was applied to previously blanched mushrooms.Mushrooms suffer a volume change during blanching due to the release of occluded air in the vegetable tissue.McArdle and Curwen [22] have shown that the shrinkage reaches up to 25%, resulting in a more compact matrix.As a consequence the thermal properties during the freezing process of previously blanched mushrooms will not take into account the air as a component of the blanched mushrooms.The specific heat, thermal conductivity and density of the mushroom (between -40 and 20°C) was considered dependent with temperature.The typical composition of the Agaricus bisporus mushrooms considered to estimate the thermal properties were: 90.5% moisture content, 4.9% carbohydrates, 0.3% fat, 3.5% protein, and 0.8% ash, as given by Berna ´s et al. [23].The moisture contents (%) of the Agaricus bisporus mushrooms used in the present work were verified experimentally by drying triplicate samples in an oven at 80°C until reaching constant weight.The models proposed by Choi and Okos [24] were implemented to estimate the thermal properties as a function of temperature and composition of the foodstuff.The thermal conductivity was: where k is the global conductivity, k i is the thermal conductivity of the component i (where i corresponds to the different components: water, ice if the temperature is lower than the initial freezing temperature T f , carbohydrate, fat, etc.), x i v corresponds to the volumetric fraction of each component.
The density of the product was calculated using: where q(T) is the global density and q i is the density of the component i (water, carbohydrates, ice, ash, etc.).The fractions ''x i '' corresponds to the mass fraction of each component.
The specific heat of the mushroom was estimated using the following Eq.17 according to the published work by [25]: The ice content as a function temperature (at T \ T f ) was estimated using the equation proposed by Miles et al. [26]: where x ice is the mass fraction of ice, x w is the mass fraction of water in the foodstuff, and T f is the initial melting point of the product.The initial freezing point of the mushroom obtained from the freezing curves using the tangent method was -1.2°C [27]; this value was in agreement with initial freezing temperatures for mushrooms values reported in [28]; T f = -1.8± 0.8. Figure 3a, b and c show the functional relationships of the thermo-physical properties with temperature.The volumetric specific heat in Eq. 17 was obtained by multiplying the q(T) by the Cp(T), and then numerically integrated using Eq. 4 to obtain the enthalpy function (Fig. 4a).The same procedure was carried out with the k(T) in order to obtain the Kirchhoff function (E(T)) according to Eq. 5 (Fig. 4b).Figure 4c shows the function E(H).This information was used as input in the numerical code in order to simulate the freezing process.
One assumption of the model used to calculate the thermal properties was to consider blanched mushrooms as an homogeneous material, without air present in the foodstuff.
During the blanching process a contraction takes place due to removed air which was occluded in the fresh mushroom structure.

Experimental procedure
The freezing experiments were applied to previously blanched whole Mushrooms, which were afterwards placed over a metal mesh in a tunnel freezer at -15°C.The average dimensions of the sliced and whole processed mushrooms were 0.0075 ± 0.005, 0.033 ± 0.004 and 0.026 ± 0.003 m in thickness (e) (sliced mushrooms), button (cap) diameter (D) and total length (L), respectively.The blanching procedure was carried out, according to Lespinard et al. [5], during 7 min at 90°C in a water bath.One group of the blanched whole samples were sliced and introduced in a tunnel (FRIOTECNOLOGI ´A, Argentina) placed over a metal mesh in order to be frozen individually at an external air temperature of -15°C.Thermocouples type T (copper-constantan (Cu-CuNi)), inserted in the mushrooms, were connected to an acquisition system (Testo 175, Testo AG, Germany) in order to record the time-temperature evolution in the product and cooling medium every 10 s for several samples at each run.
The cooling air velocity in the tunnel freezer was measured using a hot wire anemometer (TSI model 1650); the value obtained was v = 2.5 m s -1 .
Freezing of whole and sliced mushrooms were carried out in triplicate with the same operating conditions and similar mushroom sizes.In order to estimate the surface heat transfer coefficient the transient method was used; the numerical solution of heat conduction equation is the most appropriate method when dealing with heterogeneous foodstuffs, complex 3D geometries or variable boundary conditions [9,29].To estimate the heat transfer coefficient, a bronze mushroom shaped object was manufactured.Bronze was chosen as a test material due to its high thermal diffusivity, which assures an almost instantaneous uniform temperature profile.A thermocouple was inserted at the centre of object to sense the time-temperature history during cooling in the tunnel freezer, which operated under the same conditions as the mushrooms freezing experiments (previously described).The following thermo-physical properties for bronze were used: q = 8,470 kg m -3 , C p = 376.81J kg -1 K -1 , k = 122.87W m -1 K -1 .Different heat transfer coefficients were used to simulate temperature profiles; experimental and predicted temperatures for each proposed h coefficient were compared.The heat transfer coefficient that minimized the residual sum of squares (RSS) given by Eq. 19 was selected.
The h value obtained was in agreement with h values calculated with the equation proposed by Earle [30], that considers air velocities v \ 5 m/s.The h value was similar (h = 15.5 W/m 2 °C) to the value calculated in this work.
During the freezing process the rate of heat transfer from the food to the external medium depends on several factors such as the velocity of the nearby fluid, the surface roughness, the geometry of the food, etc. However these factors are very difficult to be accounted for in the mathematical formulation of the problem.In order to overcome this limitation the multiple effects are considered and quantified by the surface heat transfer coefficient (h) [31].The experiments to determine the surface heat transfer coefficient assumed an effective h value, therefore the boundary condition was considered with the same h for the entire surface exposed to the external medium; this represents a limitation to the present finite element code.A great number of publications apply an overall surface heat transfer coefficient for freezing simulations [32][33][34][35][36][37][38][39][40], however a more rigorous procedure would be to calculate the surface heat transfer coefficient by solving the Navier-Stokes equation locally.This approach has been done by Trujillo [41], Trujillo and Pham [42,43], which have used the CFD method and the FLUENT software to predict airflow and temperature fields around a beef side and therefore the heat transfer coefficient in several places throughout the carcass surface, also Harris et al. [44] measured values of h at different positions in a lamb carcass.We understand the value and importance of these measured h values implemented specially in products with much larger dimensions and more complex shapes where there is also the possibility of water evaporation from the product.A future implementation of a variable h value throughout the surface of the foodstuff coupled with the main freezing numerical model using the Enthalpy and Kirchhoff transformation is a work yet to be developed by the authors.

Results and discussion
The numerical model which consisted in a pre-processing, main, and post-processing routines, were all coded in Matlab language.The codes implemented in 2D and 3D geometries were validated using the analytical solutions of the heat conduction equation with convective boundary conditions [19,45].The accuracy and convergence of the numerical predictions were corroborated, calculating the difference between numerical and analytical solutions [46,47].The numerical results were in agreement with the analytical solutions for all systems studied (sphere and cylinder), therefore the numerical code was considered theoretically validated.
The numerical code implementing the Enthalpy and Kirchhoff transformations was compared with experimental freezing curves found in literature for irregular minced meat objects containing 75% water [33,34,36] since there is no analytical solution to the heat transfer equation when the thermal properties are temperature dependent.The numerical predictions resulted in agreement with the experimental data, even though the material exhibited a sharp phase change.Furthermore the code was compared to time-temperature predictions during the freezing process in mushrooms under industrial operating conditions.

Comparison of the numerical model with experimental results for the freezing process
The numerical simulations obtained with the code, were contrasted with experimental time-temperature freezing curves of mushrooms.Freezing experiments were undertaken in order to verify the accuracy of the solutions obtained by the finite element formulation.It can be seen from the results that a good agreement was obtained between the experimental measurements and the numerical predictions.The absolute maximum difference (infinite norm of the error = kT exp -T sim k a ) for all experiments considering all runs was less than 3.10°C.
As an example Fig. 5a and b show the predicted and experimental temperatures for the whole and sliced mushroom shape (2D axial symmetric domain).
Additionally the sliced mushroom which corresponds to a 3D geometry was analyzed considering a one-directional heat transfer flow, i.e. the heat flow (q = rTÁn) perpendicular to the z axis was considered not significant.In order to compare both results a one-dimensional phase change problem taken into account the thickness of the mushroom (e = 0.0072 m) was coded using Matlab 6.5, neglecting the x and y heat transfer contribution (see Fig. 6).The time temperature evolution was compared in z = e/2 = 0.0036 for the one-dimensional problem and for the 3D problem the internal point ((x, y, z) = (-0.001,0.0245,0.0036) in meters), which have the same z coordinate.Figure 5b shows the experimental measurements and the predicted temperature for the 1D and 3D models.As can be seen the experimental freezing curve satisfactorily agrees with the three dimensional model proposed, therefore it can be concluded that the freezing of sliced mushrooms cannot be approximated as a one-dimensional problem.The object of using a 1D approach is to demonstrate that there is a heat transfer contribution in all directions (x, y, z) and that sliced mushrooms cannot be approximated as a one dimensional object where only the thickness of the product is considered (see Fig. 6).Freezing simulations are complex problems to be solved mathematically, especially 3D problems, there is however a number of foods that can be approximated as an infinite slab that have been simulated in the past, such as hamburgers, another can be pastry discs dough, but this is not the case for mushrooms due to its average dimensions and shape.If a sliced mushroom shape is considered a one-dimensional heat flow problem, the calculated freezing times would be overestimated with respect to the actual 3D model, resulting in a much higher energy demand and quality loss.
Figure 7a shows the temperature distribution after 972 s inside the tunnel freezer and Fig. 7b the temperature distribution at the symmetry plane of the irregular shape body.It can be observed that the corners of the mushrooms reached lower temperatures much faster than the thicker central zone.This observation although trivial shows that the numerical predictions over the entire domain did not oscillate or reach unrealistic temperature values, especially at the narrower corners of the foodstuff, which usually destabilizes the numerical solution.Previous works that simulated freezing of mushrooms [12] reported unrealistic results of temperature at the boundary nodes when using the finite difference method.For example the boundary temperatures were lower than the external medium temperature; another problem was the temperature oscillation at the surface of the foodstuff.This could be attributed to the approximation method and rectangular mesh implemented in finite difference; the rectangular elements are not a sufficiently good approximation when working with curvilinear and irregular boundaries.This works showed the excellent results obtained with the finite element technique which uses triangular and tetrahedron elements that have the advantage of adjusting to any irregular shaped food.Therefore the algorithm can be considered robust and convergent, independent of the functionality of the thermal properties with temperature and the irregular domain.

Application of the numerical model to predict freezing times
Since the temperature distribution within a product varies considerably during freezing process, freezing times must be defined with respect to a position.The thermal centre is generally taken as reference, which is the location where the temperature changes most slowly.The freezing time is usually definite as the time to reach a particular temperature at the warmest point (the thermal centre) [48].For freezing, a number of final centre temperatures have been used: -5, -10 and -18°C [7].
It is known that the h value depends on several factors such as the external medium, shape and surface roughness of the foods.Typical values of h were measured and reported for various methods by several authors and were summarized in [31]: for natural air (h = 6-20 W (m 2 K) -1 ), forced convection air (h = 20-90 W (m 2 K) -1 ), air impingement (h = 80-160 W (m 2 K) -1 ), immersion in liquid nitrogen (170-425) and immersion in agitated fluid (160-1,500).In this way, simulations were carried out in order to obtain the freezing times for different external fluid temperatures and surface heat transfer coefficients values, considering an initial temperature of 10°C and two typical external fluid temperatures (-15, -20 and -30°C).The time needed to reach a value of -10°C in the warmest point (coordinates: 2D r = 0, z = 0.0139 and 3D x = 0.0023, y = 0.0254, z = 0.0036) was calculated for whole and sliced mushrooms (Table 1) considering different h values (5-1,000 W (m 2 K) -1 ) according to the before mentioned operating freezing methods.
All the numerical runs were tested for their computational speed, the maximum CPU time were less than 9 min for 3D runs and 3 min for 2D axial-symmetric model using a PC Intel(R) Core(TM) 2 6300 with a processor speed of 1.86 GHz and a RAM of 2 GB.
It is noteworthy that the computational speed of the numerical simulations of the freezing process is low, due to the simultaneous Kirchhoff and enthalpy transformation.Using this combined transformation the conductance, capacitance, and convective matrices as well as the thermal load vector in the finite element formulation remained independent of the time variable.

Conclusions
Finite element codes were developed to simulate the phase change transition during freezing of an irregular foodstuff.A combined Kirchhoff and Enthalpy formulation of the heat transfer equation was applied to deal with the highly non-linear mathematical problem and to enhance the computational speed of the numerical runs, resulting in a minimum CPU.The finite element code was found to be convergent, no oscillations of the temperature predictions were encountered in the solutions.
The numerical model was compared with experimental freezing curves of whole and sliced mushrooms, applying a two-dimensional axial-symmetric domain and 3D geometry, respectively.For both food shapes the numerical runs agreed well with experimental time-temperature curves.The freezing of sliced mushroom was also compared with a one-dimensional heat flow model, in order to verify that the mushroom slice object must be approximated as a 3D heat transfer problem.
The code was applied to obtain freezing times of mushrooms considering different operating conditions.This information can be of use for food processors to obtain valuable information to design freezing equipment and optimize processes.

Fig. 1 a
Fig. 1 a Digital photograph and characteristic dimensions of the whole mushroom and b 2D axial symmetric domain discretized into triangular elements

Fig. 2 a
Fig. 2 a Digital photograph and characteristic dimensions of the mushroom slice and b 3D mesh using tetrahedral elements

Fig. 3
Fig. 3 Thermo physical properties of the mushrooms used for the freezing process as a function of temperature a thermal conductivity b density and c specific heat Fig. 4 Functional relationships used in the combined formulation of the freezing process a Enthalpy versus temperature b Kirchhoff function versus temperature c Kirchhoff function versus enthalpy

Fig. 5 aFig. 6 Fig. 7
Fig. 5 a Experimental (opened circle) and predicted (dash line) temperatures during the freezing process for whole mushrooms (2D axial symmetric model) at an internal point of the product (point (r, z) = (0, 0.019) in meters) Ti = 22.9°C and b experimental and predicted temperature values using a 1D (at x = 0.0036 m) and 3D model at an internal point ((x, y, z) = (-0.001,0.0245,0.0036) in meters), T 0 = 6.9°C.For all experiments T ext = -15°C, h = 17 W (m 2 °C) -1 NÞdX e is the global capacitance matrix dXsðN T h r T ext ÞddX s is the global thermal load vector

Table 1
Freezing times in minutes of sliced (3D model) and whole mushrooms (2D axial symmetric model) to reach a value of -10°C in the warmest point considering different operating conditions; surface heat transfer coefficients and external fluid temperaturesh (W (m 2 °C) -1 ) T ext = -15°C T ext = -20°C T ext = -30°C