Abstract
A moving boundary model for food isothermal drying and shrinkage is applied to predict the time decay of water content and sample volume, as well as water diffusivity for chayote discoid slices in the temperature range 40–70°C. The core of the model is the shrinkage velocity , assumed equal to the water concentration gradient times a shrinkage function α representing the constitutive equation of the food material under investigation. The aim is to provide a case study to analyze and quantify differences and accuracies of two different approaches for determining the shrinkage function α from typical experimental data of moisture content vs. rescaled volume : a fully analytical approach and a shortcut numerical one.
1. Introduction
Mathematical modeling in chemical engineering historically provided the necessary support for understanding physics and transport phenomena that underlie the chemical processes. Mathematical models are also useful tools for the optimization of experimental campaigns and for scale-up from laboratory to industrial scale. During decades, more and more mathematical modeling techniques have been developed [1, 2] bridging different spatial and temporal scales from the microscopic to the macroscopic ones [2, 3] through the mesoscale [4–6]. The translation of chemical processes into equations has always been the unavoidable step for efficient plant design, process analysis, or control.
Food process engineering represents one of the most promising research fields in chemical engineering that could benefit the enhancements of the theoretical research. It was quite intentionally forgotten, in the past, due to the complexity of food materials and of related transformation processes. Among them, drying is undoubtedly one of the most investigated, probably for its huge industrial impact [7–9]. Actually, the drying process is very complex as it implies material transformation at several spatial scales (e.g., porous structure variation and volume changes).
Far from accomplishing a multiscale mathematical formulation able to account for the complexity of the all the phenomena involved [10–12], mathematical models for drying try to capture those features that are more important from the technological point of view. Shrinkage is the most relevant phenomenon connected to drying since it influences consumer quality perception, costs for transportation, and storage [13]. Among the wide literature on this topic, it is worth mentioning some works: [14] in which shrinkage kinetic laws are discussed; [15–17] in which applications of mathematical modeling techniques to food drying can be found; [10] in which a general discussion on advanced computational modeling for drying can be found; [13, 18–20] in which different advanced approaches for linking water content evolution to shrinkage can be found.
In particular, Adrover et al. [21, 22] recently developed a mathematical model with moving boundary for drying of food materials by suitably modifying a classical moving boundary model originally developed in [23, 24] for mass transport, swelling, and dissolution in polymers [25, 26]. The novelty consists in the introduction of a new constitutive equation linking the boundary movement to the water concentration gradient through a proportionality factor α representing the fingerprint of the food material under investigation.
In [21, 22], the authors were able to predict both water content evolution and volume reduction following two different approaches for the determination of the proportionality factor α from experiments: a fully analytical approach and a shortcut one. Reasonably, from a computational point of view, the first one could be more difficult to apply in many cases of practical interest, e.g., for very complex geometries of food samples, while the second one could not be accurate.
In this paper, the comparison between the two approaches is carried out by using literature experimental data [27] on chayote discoid slices. In this case, due to both the regular shape and the aspect ratio of the samples, both approaches can be easily adopted. The aim is to provide a case study to test differences in terms of accuracy of prediction of water content evolution, volume reduction, and other significant physical quantities such as the water diffusivity.
2. Moving Boundary Model for Discoid Samples
We briefly recall the model equations derived in [21, 22].
Let be a characteristic length of the food sample, ϕ the point-wise water volume fraction, the uniform initial water volume fraction, and D the water diffusivity at the operating temperature.
By introducing the dimensionless space and time variables , , , and and the dimensionless differential operators and , the moving boundary model equations for the normalized water volume fraction attain the form:where is the mass transfer Biot number, is a mass transfer coefficient, is the solid (pulp) density, is the air density at the operating temperature, is the water partition ratio between the gas and the solid phases , H being the absolute air humidity (kg water/kg dry air).
The core of the model is the shrinkage velocity , proportional to the water concentration gradient times a proportionality function tuning, at each point of the system, the relationship between water flux and volume reduction. The shrinkage velocity evaluated at the boundary controls the movement of the boundary itself. The time evolution of the sample boundary, as described by equation (3), is consistent with the classical description of boundary movement induced by the transfer of a diffusing substance across the interface [23, 24, 28]. In point of fact, equation (3) represents a generalization of a classical Stefan condition because it accounts for structural changes of the material during the drying process through the introduction of the shrinkage function .
In dealing with a discoid sample (radius , thickness ) with , we choose as a reference length so that the dimensionless initial domain is , being the discoid aspect ratio.
For high values of the aspect ratio , a one-dimensional model can be readily adopted, describing the time evolution of the rescaled water volume fraction along the axial coordinate (associated with the smallest initial dimension ) and the time evolution of the dimensionless sample thickness (uniform along the radial direction):
In this 1-d approach, both radial shrinkage and water flux from the discoid lateral surface are neglected.
3. Estimation of the Shrinkage Factor
If we adopt a 1-d model for describing the drying process of a discoid sample, the shrinkage factor can be assumed a priori or it can be estimated from the thickness calibration curve [29, 30], i.e., from experimental data of rescaled thickness vs. the moisture ratio as follows:where is the rescaled water volume fraction evaluated, at each time instant t, in a suitable point , called probe point P, placed on the sample surface and evolving in time together with the surface itself. The moisture ratio (or rescaled moisture content) is defined as , where is the amount of water at time t, is the initial amount of water in the sample, W is the sample weight (water + pulp) at time t, and is the dry sample weight. In terms of dimensionless variables the moisture ratio is ) .
In Adrover et al. [21, 22], we have shown that a proper choice of the probe point is a point exhibiting the maximum displacement (shrinkage). In this 1-d problem, since ψ is exclusively a function of ζ, the probe point is necessarily the surface point located at and evolving in time together with sample thickness .
From equation (10), it is evident that the thickness calibration curve G and, more specifically, its derivative actually furnish an experimentally derived shrinkage factor :that depends on the integral quantity and not on the required probe point concentration .
In order to recover from , it is necessary to identify a function relating to , thus obtaining
By considering that explores the entire range of values of ψ (from at to for ), once has been derived, we can adopt the same expression for to evaluate the shrinkage velocity at each point in the domain.
The simplest model that can be adopted for the function is a linear model, i.e., that implies on the boundary and consistentlyat each point in the sample domain.
By observing thatit is easy to see that the linear approximation underestimates at short time scales, when , while it overestimates at large time scales, when and (Figure 8 in Appendix). However, the linear approximation , independent of , and the resulting shrinkage function equation (13) represents an acceptable compromise between simplicity and a reasonable physical description of the drying process. Moreover, equation (13) represents a good starting point for a more accurate estimate of the shrinkage factor .
In order to derive a more accurate explicit expression for , we can adopt two different strategies: (i) a fully analytical approach [21] that can be easily applied for foods characterized by linear of quadratic calibration curves and (ii) a shortcut numerical approach [22] that can handle any nonlinear function .
In order to compare the two different approaches, we focus on experimental data of convective hot-air drying of chayote discoid samples characterized by a high initial aspect ratio so that the 1-d model, described above, can be reasonably applied.
4. Chayote Discoid Samples Air-Drying
We analyze experimental data of convective hot-air drying of chayote discoid samples (data from [27]). Fruits were washed and peeled. Cylindrical slices with initial radius mm and initial thickness mm were prepared. The discoid aspect ratio is . Chayote drying (at ) was carried out in a convective dryer with air velocity 2 (m/s). The initial moisture content was g water/100 g product. The initial water volume fraction was [31]. Available experimental data are the moisture ratio vs. time (min) and and the thickness calibration curves vs. at the four temperatures analyzed.
In [21, 22], Adrover et al. have shown that low values of lead to a smoother and flat boundary profile of the sample, while higher values lead to pronounced cusps. By considering the air velocity in the convective dryer and since no information are reported in [27] regarding significant nonuniformities of the sample thickness along the radial direction, we assume in the further analysis a low value for , namely, .
Figures 1(a) and 1(b) show experimental data for vs. and two different approximating functions, both valid in the entire range of temperature analyzed: a second-order polynomial function (Figure 1(a)) and a more accurate fourth-order polynomial function (Figure 1(b)):
(a)
(b)
The fourth-order polynomial function better captures the initial linear behaviour (for large values of ) as well as the significant increase of the volume reduction rate for larger time scales (small values). Although both functions and approximate well the global behaviour ( for and for ), the two different approximations lead to very different experimental shrinkage functions as shown in Figure 2.
We adopt the quadratic function for applying the fully analytical approach and the more accurate fourth-order function for the shortcut numerical approach. The final goal is to estimate with both approaches the effective water diffusion coefficient D and its dependence on temperature T and to compare the results.
It should be observed that there is no physical difference between the shrinkage function obtained with the fully analytical approach and that obtained from the shortcut approach. The two approaches furnish different approximations of the shrinkage function that is an intrinsic property, a sort of fingerprint, of the material under investigation.
In point of fact, during dehydration, the material exhibits structural changes. As a consequence, sample volume can be reduced by an amount that can be greater, less, or equal to the volume of released water. The introduction of the shrinkage function allows us to take this feature into account because α represents the proportionality factor between the pointwise shrinkage velocity and the local water flux. This proportionality factor changes with the water content in a nonlinear fashion for many materials as can be deduced, in a straightforward way, from the nonlinear behaviour of or, better to say, from the nonconstant behaviour of .
For a material exhibiting a linear calibration curve , the corresponding shrinkage function is constant . If , the material exhibits ideal shrinkage because volume reduction equals, at each time instant, the volume of released water.
However, in most cases, is a nonlinear function of . In the specific case of the chayote, is an increasing nonlinear function of , better approximated by a fourth-order polynomial if we want to accurately describe its behaviour for small values of , that is, on long time scales of the dehydration process. This implies that while sample volume reduction and volume of released water are comparable with each other at short time scales, on the contrary, on longer time scales, when the rescaled moisture content is below 40%, a structural collapse of the sample occurs and the volume reduction is two to three times larger than the volume of release water (Figure 2).
5. Analytical Approach
The analytical approach can be easily applied when the calibration curve is well approximated by a linear or a quadratic function because the first step is to enforce equation (14) that can be rewritten as
By solving equation (17) with respect to , we obtain an analytic expression for as a function of .
For a linear G function , equation (17) reads asthat can be solved with respect to , thus obtaining
Analogously, for a quadratic G function, e.g., equations (15) and (17) read asthat, solved with respect to , gives the following expression for :
In the case of the calibration curve being an n-th order polynomial, then the identification of the function would require the solution of an n-th order equation in and this would be extremely difficult if not analytically impossible.
The second step in the analytical approach is to relate to by adopting the following model that has been specifically derived for a one-dimensional drying and shrinkage process (Appendix) and that explicitly takes into account the influence of as follows:where
The symbol represents the smallest positive root of equation (26), and is the asymptotic rescaled thickness. The analytic derivation of the equations (22)–(26) relating to and is reported in Appendix and is independent of the shrinkage function .
By making use of the analytic expression , equations (19) or (21) (the latter in the case of chayote discoid samples), and the analytical expression for , equations (22)–(26), the more accurate expression for the shrinkage factor is obtained:
This can be used to solve the moving boundary model equations (6)–(9).
Figure 3(a) shows the behaviour of obtained with the linear approximation , i.e., (equation (13), red dashed line) and with the nonlinear model equation (27) (blue curve) for chayote discoid samples with thickness calibration curve (equation (15)).
(a)
(b)
Figure 3(b) shows the corresponding 1-d moving boundary model predictions for vs. and the comparison with experimental data. We observe that the model prediction for vs. with obtained with the linear approximation (red dotted curve) is already satisfactory (), but the improvement in model predictions obtained by adopting the nonlinear model equation (27) is significant (blue curve, ) and leads to a more accurate estimate of water diffusivity.
6. Shortcut Numerical Approach
In Section 5, we have shown that the analytical approach requires two steps in order to evaluate an accurate expression for the function, directly relating to and . The first step is to evaluate the function relating to . The second step is to evaluate the function relating to and .
The basic idea behind the shortcut numerical approach is to make a direct use of numerical data obtained with the moving boundary model in order to obtain an accurate expression directly for .
The shortcut numerical approach stems from the following observation. In Section 5, we have shown that the shrinkage function (equation (13)) obtained with the linear approximation gives reasonably good results in terms of model predictions of experimental data for vs. . Figure 4 supports this observation as it shows the good agreement between the experimental thickness calibration data for chayote discoid samples and 1-d moving boundary model predictions (red dotted curve) with given by equation (13) adopting, in this case, the fourth-order more accurate approximating function (equation (16)).
Therefore, the first step in the shortcut approach is the numerical integration of the moving boundary model with the shrinkage function equation (13) obtained with the linear approximation (and with the more accurate function) in order to obtain a numerical accurate estimate of the function, i.e., (Figure 5(a), dots).
(a)
(b)
The second step is to obtain an analytic expression for the function by a least-square best fit of numerical data vs. by means of the versatile sigmoid function:where , , and are best fit values (Figure 5(a), continuous line).
The more accurate expression for the shrinkage function reads as
Figure 5(b) shows the significant difference between the shrinkage function obtained with the linear approximation , i.e., (equation (13), red dashed line) and with the nonlinear model equations (28) and (29) (blue curve) for chayote discoid samples with thickness calibration curve (equation (16)).
Figure 4 shows the corresponding 1-d moving boundary model predictions for vs. and the comparison with experimental data. Like for the fully analytical approach, we observe that the model prediction for vs. with obtained with the linear approximation (red dotted curve) is quite satisfactory but the adoption of the nonlinear more refined model for , equations (28) and (29), permits the model to perfectly reproduce the experimental thickness calibration data in the entire range of values.
7. Estimation of Water Diffusivity
All the above analysis is independent of the water diffusivity D as the thickness calibration curve vs. is time independent.
The moving boundary model is a time-dependent model usually written in terms of dimensionless space and time variables. Specifically, the dimensionless time adopted is , being the initial sample thickness, i.e., mm.
Once the shrinkage factor has been determined, resulting the same for the four temperatures under investigation, since no significant differences have been observed in the thickness calibration curves at the different temperatures, the moving boundary model furnishes a unique time-dependent curve as well as a unique time-dependent curve vs. τ for the four temperatures.
The value of water diffusivity at each operating temperature is determined by means of a least-square best fit of the model curve onto the experimental data for time decay of the moisture ratio vs. t at the corresponding temperature T, for the same value of , namely, , adopted for the estimate of the shrinkage function .
We followed two different approaches for estimating the shrinkage function : the analytical and the shortcut numerical approaches leading to two different expressions for , namely, (equation (27)) for the fully analytical approach and (equation (29)) for the shortcut numerical approach, and therefore to two different time-dependent curves , and , respectively.
Figures 6(a) and 6(b) show the comparison between experimental data for vs. (t) (min) for chayote discoid samples and the two model curves.
(a)
(b)
and with the best fit diffusivity values are shown in Figure 7. Specifically, from the analytical approach, one obtains (m2/s), (m2/s), (m2/s), and (m2/s), with an Arrhenius correlation function with (m2/s) and (K). From the shortcut numerical approach, one obtains (m2/s), (m2/s), (m2/s), and (m2/s) with an Arrhenius correlation function with (m2/s) and the same activation energy (K).
The agreement between model curves and experimental data for the time decay of the moisture ratio is extremely good for both approaches, and differences between estimated diffusivities are below 3% for all temperatures analyzed.
Diffusivities estimated from both approaches are comparable with literature data for other vegetable products such as carrots ( (m2/s), 60–90°C), mango ( (m2/s), 40–70°C), and potatoes ( (m2/s), 40–85°C) in which shrinkage has also been considered [32–34].
This result confirms that both approaches are robust and reliable and can be indifferently adopted for predicting the time decay of the moisture ratio as well as sample shrinkage during the entire drying process.
8. Conclusions
In this paper, the comparison between the two approaches developed in [21, 22] for mathematical modeling of food drying with shrinkage is carried out.
The developed mathematical model consists of a mass balance equation for water volume fraction evolution, coupled to an equation for the movement of the boundary. The relationship between water concentration and boundary movement is introduced into the system through the dependence of a shrinkage velocity on the concentration gradient times a shrinkage factor α. This represents the constitutive equation of the material.
The two developed approaches (fully analytical and shortcut), precisely dealing with the determination of the shrinkage factor α, were analyzed here and applied to the literature data of chayote discoid sample drying.
We showed that, for the case discussed here, both approaches are able to provide accurate predictions of the experimental data for moisture fraction and thickness reduction, as well as a good estimation of the physical quantities involved (first of all, the water diffusion coefficient) with almost the same computational efforts.
However, this last consideration is true only in this case. The adoption of a one-dimensional simplified model is, in fact, due to conditions of the samples (the symmetrical shape and the value of the aspect ratio) which, as an instance, are in general not found in industrial applications.
For such complex cases, although it is reasonable to assume no differences in terms of prediction capabilities, the computational efforts connected to the fully analytical approach could be very huge compared to the shortcut one or, worst, it could not be pursued at all. Hence, the adoption of the shortcut approach, having demonstrated its efficiency here, can be an attractive opportunity (the only one in the worst cases) to have an accurate description of the evolution of the system through a mathematical model based on balance equations.
Appendix
In this appendix, equations (22)–(26) are derived. Equations (22)–(26) relate to , P being the probe point located at , .
The model equations (22)–(26) are derived from the analysis of vs. in the two limiting cases of short and long (asymptotic) time scales and are independent of the shrinkage function .
A. Short Time Scales
At the beginning of the process, when the shrinkage is small and therefore negligible, it can be assumed that the water concentration is almost uniform in the inner part of the sample and equal to the initial uniform value . Moreover, we can assume that the concentration gradient is localized in a boundary layer of thickness δ close to the boundary where the water volume fraction is .
By assuming a linear concentration profile in the boundary layer, the average water concentration in the sample reads as
The thickness δ of the boundary layer depends on and can be estimated by enforcing the linear concentration profile in the boundary layer and the boundary condition:
By substituting equation (A.2) for δ into equation (A.1) for , one obtains an analytic expression for as a function of and :that is, valid at short time scales, i.e., when .
By expanding equation (A.2) in power series about the point , one obtainsand therefore, a quadratic behaviour is predicted for vs. at the short time scales of the drying process, i.e., for .
B. Long Time Scales
When approaching the asymptotic behaviour, the sample shrinkage has almost reached its asymptotic value and concentration gradients are very small. Since the shrinkage velocity is proportional to the pointwise concentration gradient, its contribution becomes negligible and therefore the water volume fraction profile can be evaluated by solving the following pure diffusion problem in a 1-d fixed boundary domain :where attains the formand the quantities and can be evaluated as
In the asymptotic limit , given the exponential decay of , the leading term is the zero-order term associated with the dominant eigenvalue , given by the smallest positive root of equation (B.4), so that
By merging together equation (B.6), the following relation between and is obtained, valid only for long time scales:
Equation (B.7), rewritten in terms of the original variables and , reads ashighlighting a linear relationship between and at long time scales of the drying process, i.e., for .
C. A Connection between Short and Long Time Scales
In order to have an expression for vs. reliable in the whole range , the following sigmoid function is adopted (equation (22) in the main text)that satisfies the two physical constraints (uniform initial water concentration) and (uniform water concentration equal to the equilibrium concentration in the asymptotic limit).
The two parameters and entering equation (C.1) can be obtained by enforcing the two limit behaviours, equation (A.4) for and equation (B.8) for .
By expanding equation (C.1) in power series about the point , one obtains
By comparing equation (C.2) with equation (A.4) valid when , one obtains for :
By expanding equation (C.1) in power series about the point , one obtains
By comparing equation (C.4) with the asymptotic behaviour equation (B.8), one obtains the following expression for :where is given by equation (B.7).
It is important to point out that the analytical expression derived above is independent of the shrinkage factor and only the asymptotic rescaled thickness enters into equations (C.1), (C.3), and (C.5). Therefore, it can be applied, in principle, to any food material characterized by a dominant one-dimensional shrinkage.
The analytic derivation proposed above is the one-dimensional version of the two-dimensional relation derived in [21] for 2-d square cross sections.
Reliability and accuracy of the analytical expression equations (C.1), (C.3), and (C.5) containing no fitting parameters are confirmed by direct comparison with numerical results for vs. , shown in Figure 8(a), obtained in the ideal shrinkage benchmark case , , for . Figure 8(b) shows the good agreement between numerical data for vs. and the model predictions in the ideal shrinkage case for which equation (19) holds true with .
(a)
(b)
Nomenclature
| : | Biot number (−), equation (4) |
| D: | Water diffusion coefficient (m2/s) |
| : | D estimated from (m2/s) |
| : | D estimated from (m2/s) |
| : | Function relating to (−), equation (12) |
| : | Thickness calibration curve (−), equation (11) |
| : | 2nd order polynomial approx. of (−), equation (15) |
| : | 4th order polynomial approx. of (−), equation (16) |
| : | Dimensionless sample thickness (−) |
| : | Initial sample thickness (m) |
| : | Reference length (m) |
| : | Function relating to (−), equations (19) and (21) |
| : | Initial sample radius (m) |
| S: | Dimensionless sample surface (−) |
| : | Sigmoidal function (−), equation (28) |
| T: | Operating temperature (K) |
| : | Dimensionless shrinkage velocity (−), equation (4) |
| V: | Dimensionless sample volume (−). |
| : | Shrinkage function (−) |
| : | Experimental shrinkage function (−), equation (11) |
| : | Polynomial expansion coefficients (−), equations (15)–(16) |
| : | Dimensionless gradient operator (−) |
| : | Dimensionless time (−) |
| ϕ: | Water volume fraction (−) |
| : | Initial water volume fraction (−) |
| : | Rescaled water volume fraction (−) |
| : | Average value of ψ (−), equation (14) |
| : | ψ at equilibrium (−) |
| : | ψ at the probe point P (−). |
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.