Analytical solution for the coupled heat and mass transfer formulation of one-dimensional drying kinetics

Accepted Manuscript Analytical solution for the coupled heat and mass transfer formulation of onedimensional drying kinetics J. García del Valle, J. Sierra Pallares PII:

S0260-8774(18)30084-0

DOI:

10.1016/j.jfoodeng.2018.02.029

Reference:

JFOE 9184

To appear in:

Journal of Food Engineering

Received Date: 29 November 2017 Revised Date:

22 February 2018

Accepted Date: 27 February 2018

Please cite this article as: García del Valle, J., Sierra Pallares, J., Analytical solution for the coupled heat and mass transfer formulation of one-dimensional drying kinetics, Journal of Food Engineering (2018), doi: 10.1016/j.jfoodeng.2018.02.029. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

J. Garc´ıa del Vallea,∗, J. Sierra Pallaresb a Area ´

RI PT

Analytical solution for the coupled heat and mass transfer formulation of one-dimensional drying kinetics

SC

de Energ´ıas, Facultad de Ingenier´ıa Civil y Mec´ anica , Universidad T´ ecnica de Ambato, Avd. Los Chasquis y R´ıo Payamino s/n, Ambato (Ecuador) b Departamento de Ingenier´ ıa Energ´ etica y Fluidomec´ anica, Escuela de Ingenier´ıas Industriales, Universidad de Valladolid, Paseo del Cauce 59, Valladolid (Spain)

M AN U

Abstract

In this investigation the analytical solution for the coupled heat and mass transfer problem applied to the drying kinetics of thin layer slabs is presented. The coupling of the conservations laws are boundary based, while the inner matrix formulation is decoupled. The water driving mechanisms are assumed to be liquid capillarity in the inner matrix and convective evaporation on the interface. A number of physically possible boundary conditions are described, being the

D

analytical solution obtained for two commonly encountered cases in foodstuff drying experiments. Two distinct drying stages are derived as a result of the

TE

formulation of the problem. The first one is characterized for the influence of the convective coefficient as well as the ambient relative humidity. The second

EP

one is controlled by the internal effective diffusivity. In the first part of this study, the solutions have been obtained in pure theoretical grounds. The second part deals with the model validation by application to experimental data

AC C

found elsewhere.

Keywords: Drying kinetics, analytical solution, Fourier series, Chilton-Colburn analogy

∗ Corresponding

author Email addresses: [email protected] (J. Garc´ıa del Valle), [email protected] (J. Sierra Pallares)

Preprint submitted to Journal of Food Engineering

February 28, 2018

ACCEPTED MANUSCRIPT

1. Introduction

RI PT

Thin layer drying under a convective atmosphere is key process in a number of industries from construction materials to dehydrated food. Heat and mass transfer in a porous media is the underlying physical problem in these 5

applications, being the associate mechanisms described and classified in (Datta, 2007). For dehydrated foodstuff at low temperature, which is the subject of this

SC

investigation, unsteady capillarity flow is the relevant mechanism of moisture transport within the food matrix.

Exhaustive formulations of drying kinetics have been proposed in the past (Luikov, 1966) and (Whitaker, 1998). Although based on different approaches, a

M AN U

10

similar system of partial differential equations is obtained based on a number of model constants related with not easily measurable physical properties. Numerical solutions for such formulations have also been presented, (Ilic and Turner, 1989) and (Liu and Cheng, 1991). However, the reader shall consider here that 15

the majority of available drying measurements simply consist of time series for

D

the sample weight. This kind of measurement is of integral type, whereas the required data for complex model validation would be a local one, that is, the

TE

concentration variation for every point in the domain. Unfortunately this relation is rare to find. Exceptions exist, (Chiang and Petersen, 1987), where gamma ray densitometry is used to measure the water concentration inside the sample.

EP

20

On the opposite side, simplified models in which heat and mass transfer are

AC C

decoupled, based on the work of (Crank, 1975), have been used in the literature with different degrees of complexity. The general formulation involves an

25

exponential decay in which the exponential factor is assumed to be an effective diffusivity (“Def f ”). Arrhenius like temperature dependence for “Def f ” is

usually sought for, (Zhu and Shen, 2014) and (Bezerra et al., 2015). In addition, the effect of shrinkage has also been studied for Crank based formulations, (Ruiz-L´ opez et al., 2012).

30

In between this two, that is the binome complex-numerical and simplified-

2

ACCEPTED MANUSCRIPT

analytical, there is a number of models for which, although the heat and mass transfer are still coupled, some simplifications are introduced in relation with

RI PT

the former. The work of (Zhang and Datta, 2004) provides valuable insights in this regard. It is proposed that the evaporation mainly takes place at the 35

boundary, which implies that the mechanism of moisture transfer in the inner

matrix is liquid phase capillarity alone. This is a huge step towards easing the

SC

problem, since the heat and mass transfer is decoupled in the inner domain while remaining linked at the boundary. Similar arguments are found elsewhere, (Belleghem et al., 2014), where mathematical difficulties are found if vapour phase capillarity is considered for saturated materials (that is the case for drying of

M AN U

40

foodstuff, at least in the initial period). A solution for this model based on a numerical approach has been studied in the literature, (Tzempelikos et al., 2015a). The use of a numerical approach facilitates the implementation of temperature dependant transport coefficients and thermodynamic properties while 45

the physical insights are somehow blurred by the numerical implementation itself. In this regard, analytical solutions are preferred. In (Barati and Esfahani,

D

2013) a first attempt to solve the problem is presented, though the assumption of high Biot mass number effectively decouples the problem at the boundary.

50

TE

However, a particular solution was presented in (Garc´ıa-Alvarado et al., 2014). In this investigation, analytical solutions are presented for the coupled heat

EP

and mass transfer problem with a new approach related with the boundary and initial condition implementation. Details of this approach are provided prior to the mathematical solution. Furthermore, a thorough analysis of the influence

AC C

on the drying pattern of different model parameters has been performed based

55

on experimental data from elsewhere.

2. Analytical model In this section, the formulation of the problem is first discussed and later

analytically solved. The following assumptions are considered: 1. The domain is considered one dimensional. 3

ACCEPTED MANUSCRIPT

60

2. Shrinkage is neglected. 3. The effective diffusivity is constant throughout the entire process.

RI PT

4. Internal evaporation is assumed to be negligible, thus, the mass and heat transfer conservation laws may be written according to Eqs. 1 and 2

per volume. ∂C ∂2C = Def f · ∂t ∂x2   ∂T ∂2T ∂ ∂C ρ · cp =k· + Def f · · hl ∂t ∂x2 ∂x ∂x

SC

respectively. The concentration “C” is defined as the total mass of water

(1) (2)

M AN U

It shall be noted that while the mass law is purely diffusive, the energy balance takes into account the liquid water energy flux. An order of magnitude analysis between the conduction term and the energy flux through mass diffusion is presented in Eq. 3. Typical values have been used for the relevant properties. Provided the heat conduction term is predominant, Eq. 2 may be simplified into Eq. 4.

(3)

ρ · cp

(4)

TE

∂2T ∂T =k· ∂t ∂x2

D

k Conduction 1 = 1 ≈ −9 Mass diffusion D · (C0 − C∞ ) · cp 10 · 1000 · 4180

5. The boundary conditions depend on the particular experimental set-up.

EP

A number of cases are described in Table 1. Two different configurations may be found (1) For some dryers, samples are positioned in a solid metal65

lic tray and suspended inside the dryer chamber at a constant temperature

AC C

“T∞ ”. As a result, convection of both heat and mass occurs in the upper surface while, for the lower surface, there is negligible mass transfer. The heat transfer is somehow problematic. Depending on the thermal diffusivity of the tray and the sample distribution, either a constant tempera-

70

ture (case “A”), convective heat transfer (case “B”) or thermal resistance through the tray (case “C”) may be considered. (2) In other dryers, samples are suspended in perforated trays. Common permeability rates for

4

ACCEPTED MANUSCRIPT

the trays are between 50-80 %. A realistic approach is to assume convective boundary conditions both for the heat and mass transfer in either face, case “D”.

RI PT

75

Furthermore, cases “B”, “C” and “D” have two sub-variants, that is to consider different convective coefficients for the lower and upper faces. Table 1: Common types of boundary conditions found in drying experiments.

Heat

Upper face

Mass

SC

Lower face

Case

Heat

Mass

=0

Convective

Convective

T = T∞

“B”

Convective

∂C ∂x

=0

Convective

Convective

∂C ∂x

=0

Convective

Convective

Convective

Convective

M AN U

“A”

∂C ∂x

“C”

Combined resistance

“D”

Convective

Convective

In this investigation, solutions for boundary cases “A” and “D” with equal convective coefficients are considered. For the sake of discussion, boundary 80

case “A” will be considered hereon. The mathematical expressions for the

D

mass and heat transfer at the boundary are written in Eqs. 5 and 6 for

EP

∂C =0 ∂x

TE

the former, and in Eqs. 7 and 8 for the latter.

− Def f ·

∂C = hd · (C(T ) − C∞ ) ∂x

⇒x=0

(5)

⇒x=L

(6)

AC C

T = T∞ −k·

∂T = hk · (T − T∞ ) + hd · (C(T ) − C∞ ) · hlv ∂x

⇒x=0

(7)

⇒x=L

(8)

The boundary conditions deserve further attention. The following considerations apply: (a) The mass convective coefficient “hd ” can be related with the heat convective coefficient “hk ” by means of the Chilton-Colburn analogy, (Incropera and DeWitt, 1996), according to Eq. 9. The exponential 5

ACCEPTED MANUSCRIPT

term “n” is taken as 0.3. The air density “ρ” and the heat capacity

hk = ρ · cp · Le(1−n) hd

RI PT

“cp ” are evaluated at “T∞ ”. (9)

An expression for the dependence of the Lewis number with the temperature may be derived from (Hickey et al., 2000), Eq. 10. 1.32 − 1.12 · (T − 273.15) 100

(10)

SC

Le(T [K]) = 1.12 +

(b) Albeit the conditions at x = 0, so far the model resembles that presented in (Garc´ıa-Alvarado et al., 2014) and (Tzempelikos et al.,

M AN U

2015a). The concentration at the boundary (“C(T )”) in equations 6 and 8 is related with the saturation vapour concentration at the boundary temperature. The common practice in the literature (above references), is to assume that the actual vapour concentration is the saturation concentration multiplied by the water activity “αw ” at the temperature “T∞ ”. Two paradoxes arise from this assumption.

D

For the sake of clarity, the saturation concentration is assumed as a linear function of the temperature, Eq. 11, where “Cmin ” and

TE

“Cmax ” correspond to the saturation concentration at temperatures “T0 ” and “T∞ ”.This assumption will be further discussed later on.

85

Cmin − Cmax · (T − T∞ ) T0 − T∞

(11)

EP

Csat = Cmax +

Substituting Eq. 11 into Eq. 8 taking into account the water activity,

AC C

Eq. 12 is derived.

  ∂T hlv Cmin − Cmax = hk · 1 + αw · · · (T − T∞ )+ ∂x T0 − T∞ ρ · cp · Le(1−n) hlv hk · · (αw · Cmax − C∞ ) = hk · Λ · (T − T∞ ) + Υ ρ · cp · Le(1−n) −k·

(12)

An idea of the shape of the solution for Eq. 4, if boundary conditions 7 and 12 apply, can be easily derived from superposition. The term 6

ACCEPTED MANUSCRIPT

“Υ” would produce a linear temperature gradient from 0 at “x=0” and a function of “Υ” at “x=L”. The homogeneous solution and the

RI PT

term “hk · Λ · (T − T∞ )”, where “T” depends on the former solution,

would make the sample temperature transition from uniform temperature “T0 ” to some equilibrium temperature (wet bulb temperature),

rather than “T∞ ”, that is what is expected at the end of the drying

SC

period. This temperature drop is fueled by continuous liquid evapo-

ration due to the mass imbalance found in the boundary condition 6, rewritten as Eq. 13. This behaviour will ultimately produce negative

− Def f ·

M AN U

water concentrations.

∂C = hd · (αw (T∞ ) · Csat (TL ) − C∞ ) 6= 0 ∂x

(13)

(c) There is a second concern related with the use of the water activity in this context. While this parameter is actually the equilibrium water concentration at a given temperature, Eq. 12 uses it for unsteady computations. A similar though was briefly discussed

D

in (Datta, 2007). The physical behaviour is that the departure from the saturation concentration shall be higher at lower inner matrix wa-

TE

ter concentrations. It is argued here to implement the water activity concept in a whole different way. The idea is to assume a saturated

EP

concentration in the boundary up to the time when the concentration in the inner domain at “x=L” attains a value of αw ·Csat (T∞ ). Afterwards the problem is fully decoupled, being the boundary conditions

AC C

for the mass transfer according to Eqs. 14 and 15. ∂C =0 ∂x

⇒x=0

(14)

C = αw · Csat (T∞ )

⇒x=L

(15)

It is worth noting that this is actually the Crank solution (Crank, 1975). In this way, both the equilibrium concentration in the sample and the temperature T∞ will be reached at the end of the drying

90

process. 7

ACCEPTED MANUSCRIPT

6. Regarding the initial conditions, both the concentration and the temperature may be supposed to be uniform at “t=0”, namely “C0 ” and “T0 ”

RI PT

respectively. However, its numerical value, in conjunction with that of the surrounding atmosphere may produce different drying patterns. This

fact has passed widely unnoticed in the literature. For the sake of clarity,

95

Fig. 1 has been included in this discussion. It represents the evolution of

SC

the temperature, the saturated concentration and the inner concentration at the boundary (dotted line). Both the dryer temperature “T∞ ” and its

dew point “Tdew ” have been explicitly drawn. Furthermore, the dryer at-

M AN U

mosphere concentration “C∞ ”, the saturation concentration at the dryer

100

temperature “CT∞ ”, the sample equilibrium concentration (related with the water activity discussed earlier) “Ceq ” and the initial sample concentration “C0 ” have been considered. C

T

C0

T

a

Ceq

b

Tdew

D

CT C

nd

2 stage st

1 stage

time

TE

st

1 stage

rd

3 stage nd

2 stage

Inner concentration at x=L

b

a

a Saturation concentration at x=L

b st

nd

1 stage

2 stage st

1 stage

rd

3 stage

time

nd

2 stage

EP

Figure 1: Temperature and concentration profiles at the boundary as a function of time for cases “a” and “b”.

AC C

The following cases may be found, 105

(a) Case “a”: the initial sample temperature is higher than the dryer dew point. Water evaporation starts right from the beginning. Wet bulb temperature is reached at some point, thus keeping the surface temperature approximately constant. Water evaporation continues until the boundary concentration reaches “Ceq ” (1st stage). At this

110

point the problem becomes diffusion controlled (2nd stage) as discussed earlier, Eq. 15.

8

ACCEPTED MANUSCRIPT

(b) Case “b”: in this case the initial sample temperature is lower than the dryer dew point. Water would come into the sample, but for

115

RI PT

fresh products this behaviour is non-physical, since they are already

saturated. Thus, an initial period in which mass transfer is omitted is introduced until the sample reaches “Tdew ” (1st stage). Once here, the description provided for case “a” applies.

(c) Cases “c” and “d” are introduced as the corresponding cases for “a”

120

SC

and “b” respectively if the boundary case is changed from “A” to “D” according to Table 1.

M AN U

2.1. Analytical solution

Provided the mathematical derivation is long, it has been preferred to summarize the results in tabular form for cases “a”, “b”, “c” and “d” in Tables 2 to 5 respectively. The complete solution as well as the computer implementa125

tion (detailed in the following section) have been made available online (Garc´ıa, 2018). The following comments apply regarding the aforementioned tables:

D

1. The solution for both the temperature and the concentration are computed for each stage according the description provided for Figure 1.

130

TE

2. The relevant coefficients which appear in Tables 2 to 5 have been included in Annex A in Tables 8 to 11 respectively.

EP

3. The function “T(t)” in Eqs. 17 and 23 corresponds to the temperature solution evaluated at “x=L”, Eqs. 16 and 22 respectively. 4. The temperature solution has been omitted in the last stage due to the

AC C

fact that the heat and mass transfer is decoupled.

135

5. For cases “a” and ”b” the domain spans for whole the sample length. However, for cases “c” and “d”, since the boundary conditions are symmetrical, the domain covers half the sample.

6. The transition time “t∗ ” between stages 1 and 2 for cases “a” and “c” and stages 2 and 3 for cases “b” and “d” is obtained from “C(L, t∗ ) = Ceq ”.

140

7. The transition time “tdew ” between stages 1 and 2 for cases “b” and “d” is obtained after solving “T (L, tdew ) = Tdew ”. 9

ACCEPTED MANUSCRIPT

    T = T∞ ⇒ x = 0    ∂T BC & IC = −k · = hk · (T − T∞ ) + hd · (C(T ) − C∞ ) · hlv ⇒ x = L  ∂x     T = T0 ⇒ t = 0

T =

∞ X i=0

 2

Ji · e(−α·ϑi ·t) + Ki +

j=0

 2

Kij · e(−α·ζj ·t)  · sin(ϑi · x)

Υ x2 hk · Λ x · (L − x) · · VL (t) + T∞ − · k L k 2·L   ∞ X 2 2 1 Υ · (−1)i · e(−α·ζi ·t) − · sin(ζi · x) + 2 k L · ζi ζi · L i=0

SC

  Z x2 L 1 C = − (γ + β · T (t)) · + β · T (t) · − γ · t + β · T (t) · dt · 2 · L · Def f 6 · Def f L " ∞ ∞ X X 2 2 2 +K + Ej · e(−Def f ·ξj ·t) + Fji · e(−α·ζi ·t) + Gji · e(−α·ϑi ·t) j=1

∞ X

t∗

Z 0

−Def f ·



#

· cos(ξj · x)

Z t∗ ∂C dt = γ · t + β · T (t) · dt ∂x x=L 0

(17)

(18)

Not computed

D

 ∂C    =0⇒x=0    ∂x BC & IC = C = Ceq ⇒ x = L       C = C1st stage (t∗ ) ⇒ t = 0

C = Ceq +

∞ X

2

Oi · e(−ζi ·Def f ·t) · cos (ζi · x)

(19)

TE

Concentration

stage (decoupled)

T

mwater = A

i=0

2

Hjik · e(−α·ζk ·t)

k=0

2

(16)

 ∂C    =0⇒x=0    ∂x BC & IC = −Def f · ∂C = hd · (C(T ) − C∞ ) ⇒ x = L  ∂x     C = C0 ⇒ t = 0

M AN U

Case “a”

Concentration

1st stage (coupled)

+

+

nd

∞ X

RI PT

Temperature

Table 2: Summary of the analytical solution for case “a”.

i=0

∞  X 2 mwater Oi  = · 1 − e(−ζi ·Def f ·t) · (−1)i A ζ i i=0

EP

(20)

Temperature

    T = T∞ ⇒ x = 0    ∂T BC & IC = −k · = hk · (T − T∞ ) ⇒ x = L  ∂x     T = T0 ⇒ t = 0 ∞ X

2

C

C = C0

T

Mj · e(−α·υj ·t) · sin(υj · x)

Same as case “a”, 1st stage of drying.

C

j=0

Same as case “a”, 1st stage of drying.

T

T = T∞ + (T0 − T∞ ) ·

Not computed

C

3rd st. 2nd st.

Case “b”

AC C

1st stage (decoupled)

Table 3: Summary of the analytical solution for case “b”.

Same as case “a”, 2nd stage of drying.

10

(21)

ACCEPTED MANUSCRIPT

Table 4: Summary of the analytical solution for case “c”.

∞ X i=0

 2

Ji · e(−α·ϑi ·t) + Ki +

∞ X j=0

 2

Kij · e(−α·ξj ·t)  · cos(ϑi · x) + T∞

M AN U

 ∂C    =0⇒x=0    ∂x BC & IC = −Def f · ∂C = hd · (C(T ) − C∞ ) ⇒ x = L  ∂x     C = C0 ⇒ t = 0

(22)

SC

Υ x2 − L2 · k 2·L ∞ X
2 2 Υ · · e(−α·ξi ·t) (−1)j · cos(ξj · x) − 1 + k L · ξj2 j=1

Concentration

RI PT

T =

−

1st stage (coupled)

Case “c”

Temperature

 ∂T    =0⇒x=0    ∂x BC & IC = −k · ∂T = hk · (T − T∞ ) + hd · (C(T ) − C∞ ) · hlv ⇒ x = L  ∂x     T = T0 ⇒ t = 0

  Z L 1 x2 C = − (γ + β · T (t)) · + β · T (t) · − γ · t + β · T (t) · dt · + K 2 · L · Def f 6 · Def f L " !# ∞ ∞ ∞ X X X 2 2 2 + Ej · e(−Def f ·ξj ·t) + Gji · e(−α·ϑi ·t) + Hjik · e(−α·ζk ·t) · cos(ξj · x) j=1

i=0

k=0

BC & IC =

 ∂C    =0⇒x=0    ∂x

C = Ceq ⇒ x = L     C = C1st stage (t∗ ) ⇒ t = 0

D

2

nd

Concentration

stage (decoupled)

T

(23)

Not computed

C = Ceq +

∞ X

2

Oi · e(−ζi ·Def f ·t) · cos (ζi · x)

(24)

TE

i=0

Temperature

 ∂T    =0⇒x=0    ∂x BC & IC = −k · ∂T = hk · (T − T∞ ) ⇒ x = L  ∂x     T = T0 ⇒ t = 0 ∞ X

2

C

C = C0

T

Mj · e(−α·υj ·t) · cos(υj · x)

Same as case “c”, 1st stage of drying.

C

j=0

Same as case “c”, 1st stage of drying.

T

T = T∞ + (T0 − T∞ ) ·

Not computed

C

3rd st. 2nd st.

AC C

Case “d”

EP

1st stage (decoupled)

Table 5: Summary of the analytical solution for case “d”.

Same as case “c”, 2nd stage of drying.

(25)

Based on the solution for the concentration, the mass of water out of the

sample per surface area (“mwater /A”) is computed according Eqn. 16. This expression is equal for the four cases, albeit the particular coefficients change.

11

ACCEPTED MANUSCRIPT

145

3. Material and methods

RI PT

Details on the model implementation need further analysis before its application to experimental data. Therefore, a short account of the experimental

database and the influence of three model prerequisites are detailed in this section, namely: (1) an account of the deviation introduced by the linearization of 150

the saturation concentration, Eq. 11, (2) some considerations regarding the heat

SC

convection coefficient, (3) the influence of the number of terms of the Fourier expansion into the numerical results.

M AN U

3.1. Experimental measurements considered from the literature

Literature concerning drying experiments of sliced foodstuff is ample, however, investigations in which the air velocity, relative humidity and equilibrium concentration are reported are but a tiny fraction. The following experiments, which meet the aforementioned criteria, are considered: (Johnson et al., 1998), (Johnson and Brennan, 2000) for plantain, (Tzempelikos et al., 2015b) for peach and (Kaya et al., 2007) for quince. The boundary conditions for plantain drying

D

correspond to either model “a” or “b”, whereas the last two correspond to either

TE

model “c” or “d”. Relevant cases have been summarized in Table 6. It shall be noted that the original data is in terms of the dry based moisture content (“M”), Eq. 26, whereas the model computes the cumulative mass of water out

EP

of the sample over time per surface area, (“mwater /A”), as in Eq. 18 and 20. Furthermore, the model uses the initial and equilibrium concentration, but the initial and equilibrium moisture content is reported in the literature. In order

AC C

to uniformize the results from all the experiments, the ratio (“m/m0 ”) has been

sought for in this study instead of the moisture content. The advantage of this ratio is that the data is bound between [0,1]. Furthermore, comparison between drying curves of samples with different initial water content is eased. The relation between “m/m0 ”, “M”, “mwater /A”, “C0 ” and “Ceq ” is shown in Eqs. 27 to 28 provided that the initial density is known. M=

mwater mdry

(26)

12

ACCEPTED MANUSCRIPT

(27) (28)

RI PT

 m 1 m 1+M water · = =1− m0 1 + M0 A ρ0 · L M0 C 0 = ρ0 · 1 + M0 Meq Ceq = ρ0 · 1 + M0

(29)

Thermodynamic and transport properties have been included in Table 6. The initial density (“ρ0 ”), conductivity (“k”) and specific heat capacity (“cp ”)

SC

have been computed according the procedure proposed by (Choi and Okos, 1986), being the required coefficients taken from (ASHRAE, 2010). It must be noted that for the conductivity and heat capacity, the water mass fraction,

M AN U

(“xwater ”), used in the aforementioned model, has been taken as the average value, as shown in Eq. 30. The remaining constituents have been scaled accordingly. The initial moisture content is M0 = 1.62 for plantain, M0 = 4.25 for peach and M0 = 4.71 for quince. The equilibrium moisture content is included in Table 6. M0 + Meq 2 + M0 + Meq

(30)

D

xwater =

TE

3.2. Uncertainty introduced by the linearization of the saturation concentration The effect of this hypothesis can not be evaluated straightforward. In practice, the initial and the dryer temperature are chosen as the extremes. In a

EP

generic form, a “Thigh ” and a “Tlow ” are selected. An estimated overall error may be computed as the average absolute error between the exact concentration

AC C

“Csat ” and the linearised concentration, “Clin,sat ”, as shown in Eq. 31. The exact concentration is computed by means of the thermodynamic properties library CoolProp, (Bell et al., 2014). This is not at all correct, since the sample will spend most of the time at the wet bulb temperature, thus, a weighed average would be more precise. In any case, Eq. 31 would provide an order of magnitude of the expected deviations. Error =

    1 Pn i i C T + − T ) − C T + − T ) · (T · (T sat low low lin,sat low low high high n + 1 i=0 n n

13

(31)

ACCEPTED MANUSCRIPT

Table 6: Summary of the experimental measurements. The units for the sample conductivity (“k”) and specific heat capacity (“cp”) are [W m−1 K−1 ] and [kJ kg−1 K−1 ] respectively. T [◦C]

RH [%]

V [m s−1 ]

L [mm]

Meq

ρ0 [kg m−3 ]

k

cp

M0

40

21.4

3.6

5

0.05

1106.79

0.423

3.062

M1

50

12.8

3.6

5

0.025

1106.79

0.432

3.062

M2

60

7.9

3.6

5

0.016

1106.79

0.441

3.068

M3

70

5.1

3.6

5

0.01

1106.79

0.449

3.076

Peach

P0

40

10

2

12

0.106

1048.14

0.509

3.565

(Tzempelikos et al.,

P1

50

10

2

12

0.093

1048.14

0.520

3.570

2015b)

P2

60

10

2

12

0.080

1048.14

0.529

3.577

Q0

35

40

0.2

4

0.092

1065.07

0.513

3.611

Q1

35

40

0.4

4

0.092

1065.07

0.513

3.611

Q2

35

40

0.6

4

0.092

1065.07

0.513

3.611

Q3

45

40

0.2

4

0.084

1065.07

0.524

3.617

Q4

45

40

0.4

4

0.084

1065.07

0.524

3.617

Q5

45

40

0.6

4

0.084

1065.07

0.524

3.617

Q6

55

40

0.2

4

0.062

1065.07

0.534

3.622

Q7

55

40

0.4

4

0.062

1065.07

0.534

3.622

Q8

55

40

0.6

4

0.062

1065.07

0.534

3.622

Q9

35

55

0.2

4

0.1

1065.07

0.513

3.612

Q10

35

70

0.2

4

0.17

1065.07

0.514

3.618

1998)

Quince (Kaya et al., 2007)

155

SC

(Johnson et al.,

M AN U

Plantain (Musa)

RI PT

Name

Matrix

In order to visualize Eq. 31, Figure 2(a) has been introduced. It is a three-

D

dimensional map plot in which the color grading represent the percentage deviation according to Eq. 31. The number of terms “n” considered has been set

TE

to 50, further increments do not change the results significantly. As expected, the error increases as the temperature difference increases. Errors up to a 100 % are found for temperature differences of 60 ◦C. As mentioned earlier, pure

EP

160

numerical simulations of the whole model would be needed in order to translate

AC C

this error to the diffusivity or other model parameter. 3.3. Heat convective coefficient The uncertainties arising from the estimation of the heat convective coeffi-

165

cient are somehow difficult to evaluate. Not only the air velocity but the sample size, geometry and sample position within the tray influence its value. For case “a” and “b” a flow parallel to a flat plate may be considered, for which correlations are readily available. Moody-like flow regime diagrams for flat plates show the laminar regime extending up to Reynolds number of 1 × 106 , 14

ACCEPTED MANUSCRIPT

(White, 1995). Maximum sample holder dimensions up to 0.5 m and velocities below 5 m s−1 yield Reynolds number in the laminar range. Correlations for the

RI PT

Nusselt number in the laminar range, from which the heat convective coefficient is derived, may be found in any heat transfer textbook, (Incropera and DeWitt, 1996), as shown in Eq. 32. 1

2 · 0.3387 · Re0.5 · P r 3 N¯uf lat plate laminar =   2 !0.25 0.0468 3 1+ Pr

SC

(32)

Three parameters influence the convective coefficient, that is the sample length,

M AN U

air velocity and temperature. The latter has been proved to barely influence the results in the range 0 ◦C to 100 ◦C. Again, it has been found useful the use of a 170

three-dimensional map plot as shown in Fig. 2(b). The color grading represents the heat convective coefficient in terms of the air velocity and sample length. Air density, viscosity, Prandtl number and conductivity are taken at standard atmospheric pressure from the library CoolProp (Bell et al., 2014). It can be seen that both air velocity and sample length are determinant factors for the value of the convective coefficient.

TE

120 80 60 40 20 0 0

EP

Thigh [◦C]

100

20

40

60

80

100

450 400 350 300 250 200 150 100 50 0

Velocity [m s−1 ]

D

175

120

Tlow [◦C]

5

70 60 50 40 30 20 10 0

4 3 2 1 0 20

40

60

80

100

120

140

AC C

Length [mm]

(a) Density linearization error as a function (b) Study of the variation of “h” for a flat of the minimum and maximum temper-

plate as a function of the air velocity and

ature. The color gradient represents the

sample length at 25 ◦C. The color gradi-

mean absolute error in percentage accord-

ent represents h (W K−1 m2 ).

ing to Eq. 31

Figure 2: Study of the influence of some model prerequisites.

For case “c” and “d” a perpendicular flow over a disc shall be considered, 15

ACCEPTED MANUSCRIPT

for which fluid flow is far more complicated due to boundary layer separation

3.4. Influence of the number of terms for the Fourier series 180

RI PT

and will not be pursued in this investigation.

A C++ program has been written for evaluating the model performance. The following remarks are needed:

SC

1. Both the model and the experimental data has been transformed to the m ” by means of Eq. 27. The terms “mexp,i ” and “mmodel,i ” are ratio “ m 0

introduced for this ratio for the former and latter datasets respectively.

M AN U

2. A goodness of fit defined as the mean absolute error (“MABE”) is introduce for a measured series according Eq. 33. n

MABE =

185

1X |mexp,i − mmodel,i | n i=1

(33)

3. Transcendental equations 41 and 93 and the values of “t∗” and “tdew ” are solved by a bisection algorithm up to a relative error of 1 × 10−9 .

D

The same number of terms have been considered for all the Fourier series in the model. Initially, case “Q5 ” from Table 6 has been studied with 4 up to 200

190

TE

terms to determine the influence of the number of terms of the Fourier series into the mean absolute error, Eq. 33. Provided the procedure by which the model is

EP

fitted to the experimental data is detailed later, the best fitted model parameters have been used here and will be detailed below in this section. Results are shown in Fig. 3. It can be seen that the variation of the error is quite small indeed

AC C

for the whole range explored. From 30 terms onwards the error variation is

195

minimal, thus, 30 terms will be used in the remaining computations.

4. Results

In this section the model has been applied to the experimental data described

in Table 6 in two ways: 1. Non-linear regression in order to estimate both “Def f ” and “hk ”. 16

1

10

Nº terms

SC

RI PT

2.80 2.75 2.70 2.65 2.60 2.55 2.50 2.45 2.40 2.35 2.30

100

M AN U

Mean Absolute Error [ %]

ACCEPTED MANUSCRIPT

Figure 3: Influence of the number of terms of the Fourier series for case “Q5 ”.

200

2. Multi-parametrical influence on the drying curves of three model parameters: the convective coefficient, the dryer relative humidity and the effective diffusion coefficient. This study is a step towards understanding the

D

qualitative and quantitative influence of the aforementioned parameters,

205

TE

rather than solely influence of “Def f ” as is the case in Crank-like models. 4.1. Fitting of the model to experimental data

EP

The model is dependant on three parameters, the convective coefficient, the dryer relative humidity and the effective diffusivity coefficient. The relative humidity is the only parameter which is directly measured from the experiments.

AC C

The convective coefficient, as discussed earlier, may be inferred for existing 210

correlations. However, it has been preferred to fit both “Def f ” and “hk ” simul-

taneously to the experimental data and analyze the results. A two parameter fitting algorithm from (GSL, 2017) has been used for this purpose. The results have been included in Table 7. In addition to the effective diffusivity and convective coefficient, a number of intermediate model results have been included in

215

this table. The Biot number computed from the fitted convective coefficient has been included and will be used later in the discussion. The standard deviation 17

ACCEPTED MANUSCRIPT

(σ) for the parameters “Def f ” and “hk ” has been estimated during the fitting

of 2.7% and 3.8% respectively.

RI PT

procedure, yielding the quotient “σDef f /Def f ” and “σhk /hk ” an average value

Table 7: Numerical results for the fitting of the model to experimental measurements. The units for Def f are [m2 s−1 ], for hk are [W m−2 K−1 ] and for Ceq are [kg m−3 ]. m Def f hk MABE·102 Tdew [◦C] tdew [s] t∗ [h] m Name Case ∗ 0

Biot

“a”

9.22E-10

8.0

1.63

12.4

0.0

2.25

0.29

21.1

0.09

M1

“a”

1.15E-09

4.7

2.10

11.9

0.0

2.25

0.32

10.6

0.05

M2

“a”

1.69E-09

3.7

2.34

11.4

0.0

1.85

0.35

6.8

0.04

M3

“a”

2.69E-09

3.1

1.72

11.0

0.0

1.56

0.39

4.2

0.04

P0

“c”

4.91E-10

19.3

4.93

0.8

0.0

2.87

0.26

21.2

0.46

P1

“c”

6.81E-10

18.9

4.71

8.0

0.0

2.43

0.28

18.6

0.44

P2

“c”

1.18E-09

22.0

5.17

15.2

0.0

1.97

0.34

16.0

0.50

Q0

“c”

1.28E-10

29.5

2.06

18.5

0.0

1.09

0.25

17.2

0.23

Q1

“c”

1.40E-10

35.5

1.90

18.5

0.0

0.82

0.23

17.2

0.28

Q2

“c”

1.53E-10

44.5

1.99

18.5

0.0

0.57

0.20

17.2

0.35

Q3

“d”

2.17E-10

48.4

2.21

27.3

47.6

0.39

0.19

15.7

0.37

Q4

“d”

2.57E-10

68.8

2.24

27.3

31.3

0.23

0.16

15.7

0.53

Q5

“d”

3.12E-10

98.8

2.43

27.3

19.3

0.13

0.13

15.7

0.75

Q6

“d”

3.56E-10

55.9

2.70

36.0

79.7

0.29

0.21

11.6

0.42

Q7

“d”

4.36E-10

69.3

2.53

36.0

63.4

0.23

0.21

11.6

0.52

Q8

“d”

5.39E-10

86.9

2.40

36.0

49.5

0.18

0.20

11.6

0.65

Q9

“d”

9.50E-11

24.2

2.12

24.0

91.0

2.15

0.30

18.7

0.19

Q10

“d”

3.73E-11

14.9

1.63

28.3

416.6

4.93

0.28

31.7

0.12

TE

D

M AN U

M0

EP

220

Ceq

SC

t

The following comments may be drawn from Table 7:

AC C

1. Different model cases, based on the dew point temperature, appear in the simulation. It shall be noted here that the initial temperature is not provided in the original references, thus it was assumed to be 20 ◦C.

If different initial temperatures were considered, model cases may have

225

differed from those shown in Table 7. However, the influence on the model results of small variations in the initial temperature is largely unnoticeable.

2. The fitting error is small, with an average value for the mean absolute error defined in Eq. 33 of 2.6 %. 18

ACCEPTED MANUSCRIPT

3. The column “tdew ” represents the time at which the boundary reaches the 230

dew point temperature. Only for case “Q10 ” its value is significant. This

RI PT

feature is translated as a constant moisture content region close to the origin, as can be seen in the original manuscript, (Kaya et al., 2007).

4. The behaviour of “Def f ” is best understood if cases with the same velocity and relative humidity are considered, such as “P0 − P1 − P2 ”, “Q0 − Q3 −

Q6 ”, “Q1 − Q4 − Q7 ”, “Q2 − Q5 − Q8 ”. The effective diffusivity positively

SC

235

correlates with the temperature, as is observed elsewhere, (Zhu and Shen, 2014).

M AN U

5. Convective coefficients for cases “M 0” to “M 3” are lower than those predicted for the flat plate, Fig. 2(b), with a characteristic length of 31 mm 240

and a velocity of 3.6 m s−1 , yields a convective coefficient of approximately 41 W m−2 K−1 . An explanation may be derived if it is considered that the samples are place inside a small cup of the same diameter and depth as the sample. Due to the sample shrinkage during drying, the cup rim would likely protrude above the sample, thus reducing the velocity over it. In this situation is difficult to derive any relation for the Nusselt number.

D

245

TE

6. The influence of the velocity on the convective coefficient may be derived from cases “Q0 − Q1 − Q2 ”, “Q3 − Q4 − Q5 ”, “Q6 − Q7 − Q8 ”. It can be observed that the convection coefficient grows with the velocity. However,

250

EP

convection coefficients for cases “Qi ” are higher than those for “Pi ”, while air velocity is considerably lower in the former. This behaviour is non-

AC C

physical and requires further study, which is addressed below. 4.1.1. Analysis of the parameters of influence In order to get an insight into the model behaviour, three-dimensional color

grading maps in which the influence of two individual parameters are studied

255

independently are presented in Figs. 4(a) to 4(c) for [Def f , RH] and cases “M0 ”, “P0 ” and “Q5 ” respectively and Fig. 4(d) for [HR, hk ] for case “M0 ”. The following comments may be extracted from these figures. 1. For case “M0 ” the problem is dominated by the convective coefficient, Fig 19

ACCEPTED MANUSCRIPT

4(a), while for case “Q5 ” the problem is dominated by , Fig 4(c). Case 260

“P0 ” is equally influenced by “hk ” and “Def f ”, Fig 4(b).

RI PT

2. For case “Q5 ”, Crank based models, that is an exponential decay, would yield appropriate modelling. Since the minimization problem for “hk ”

is not well posed, the computed value is likely to differ from the actual convective coefficient, as noted earlier.

3. The Biot number from Table 7 is included in Figures 4(a) to 4(c). These

SC

265

plots are actually a way of visualizing the influence of the Biot number in the drying curve, namely (1) cases with Biot number lower than 0.1 are

M AN U

convective controlled, (2) when Biot number equals 0.5, convection and diffusion are of the same order and (3) for Biot number higher than 0.75, 270

diffusion is the predominant mechanism. Mean absolute error plots may be used as a substitute for Biot number assessment on the influence of the convection coefficient on the drying kinetics curve, (Bezerra et al., 2015). 4. From Fig. 4(d), the combine influence of “hk ” and “RH” shows opposite trends. If the “RH” increases, the drying curve flattens and if “hk ” increases, the drying curve steepen. That would not be a problem if the

D

275

TE

drying curve shape produced by each parameter were different. However, after examination of Fig. 4(d), it can be seen that instead of a point-minimum, a line-minimum does exist. This means mirror behavior

280

EP

between “RH” and “hk ”. Under these circumstances is difficult indeed to search for a minimum if both parameters are considered in the minimization algorithm. Thus, it makes more sense to either fit for [Def f , RH] or

AC C

[Def f , hk ], as is the case in the present investigation.

Further insights into the model behaviour for a particular case may be

gained if small deviations from the best fitted curve are considered for the three

285

parameters of influence. Thus, the six following combinations are explored, [2 · Def f , hk , RH], [0.5 · Def f , hk , RH], [Def f , 2 · hk , RH], [Def f , 0.5 · hk , RH] and [Def f , hk , RH ± 0.1]. The results for case “M0 ” are plotted in Fig. 5, where

combinations 1 and 2 are referred as “Def f +” and “Def f −”, respectively and 20

14

16

28

12 10 8 6 4

14

24

12 20

RI PT

12 11 10 9 8 7 6 5 4

h [kW m−2 K−1]

h [kW m−2 K−1]

ACCEPTED MANUSCRIPT

10

16

8 6

12

2 10

15

20 25 30 35 D [1 × 1010 m2 s−1]

4

40

2

3

4 5 6 7 8 D [1 × 1010 m2 s−1]

9

10

as a function of “D” and “hk ” (RH =

as a function of “D” and “hk ” (RH =

21.4 %). Nominal Bi = 0.09.

10 %). Nominal Bi = 0.46. 16 14 12 10 8 6 4 2

100 80 60 40 2.0

2.5

3.0 3.5 4.0 4.5 D [1 × 1010 m2 s−1]

h [W m−2 K−1]

120

12 11 10 9 8 7 6 5 4

M AN U

140 h [kW m−2 K−1]

SC

(a) Model mean absolute error for case “M0 ” (b) Model mean absolute error for case “P0 ”

5.0

10

15

20

25 30 RH [ %]

35

16 14 12 10 8 6 4 2 40

(c) Model mean absolute error for case “Q5 ” (d) Model mean absolute error for case “M0 ” as a function of “RH” and “hk ” (Def f =

40 %). Nominal Bi = 0.75.

9.22 × 10−10 m s−2 ).

TE

D

as a function of “D” and “hk ” (RH =

Figure 4: Model behaviour analysis. In all cases the color grading corresponds to the mean absolute error in percentage value as defined in Eq. 33.

290

EP

so on for “hk +”, “hk −”, “RH+” and “RH−”. Furthermore, the drying stage

transition at time “t∗ ” has been explicitly marked with a filled dot. The follow-

AC C

ing remarks may be drawn from this figure: 1. The effective diffusivity does not change the slope of the drying curve up to the time “t∗ ”, as can be inferred from Eq. 18 and the curves “Def f +” and “Def f −” from Fig. 5. This stage is mostly convective controlled

295

and corresponds to the first stage of drying for model “a” and “c” or the second stage of drying for model “b” and “d”. If the diffusivity increases, a higher rate of water flow to the boundary prevents the concentration from reaching “Ceq ”, thus the convective stage extends longer. In the 21

ACCEPTED MANUSCRIPT

contrary, if the diffusivity is reduced, water flow is restricted, thus quickly 300

reaching the concentration “Ceq ” at the boundary. Once the system enters

RI PT

the diffusivity controlled stage, the water removal rate is greatly reduced. 2. As expected, an increment in the convective coefficient “hk ” steepen the

drying curve at the same time that the transition time “t∗ ” is reduced. The dryer relative humidity works exactly in the opposite direction. 1.00

SC

Experimental Fitted Def f + Def f − hk + hk − RH+ RH−

M AN U

0.90

m m0

0.80 0.70 0.60

0

1

2

TE

0.40

D

0.50

3

4

5

6

7

Time [hour]

Figure 5: Sample drying kinetics variation in terms of Def f , hk and RH for case “M0 ”. The

305

EP

filled dots represent the drying stage transition at time “t∗ ”.

3. It shall be noted that up to the time “t∗ ” the drying curve is approxi-

AC C

mately a line. The explanation for this behaviour lays in the temperature behaviour. For that purpose, the sample temperature profiles, Fig. 6(a), have been included. Concentration profiles have been included in Fig. 6(b) for comparison purposes. Provided the time increment between each line

310

is 5 s, it means that within 10 min the boundary temperature has reached approximately the wet bulb temperature, which will remain constant until the time “t∗ ”. Further times are not modelled since the heat and mass transfer is decoupled. If the temperature remains constant, according the 22

ACCEPTED MANUSCRIPT

Eq. 18, the mass flow rate through the boundary also remains constant. From the concentration profile curves, Fig. 6(b), is readily visible the

315

RI PT

slowing rate of mass diffusion beyond “t∗ ” due to the piling up of the concentration profiles. 700

32 28 24 20 16 0.0

1.0

2.0 3.0 Distance [mm]

600 500 400 300 200 100

SC

Concentration [kg m−3]

36

M AN U

Temperature [◦C]

40

4.0

5.0

0 0.0

1.0

2.0 3.0 Distance [mm]

4.0

5.0

(a) Temperature profiles for case “M0 ” based (b) Concentration profiles for case “M0 ” on a time increment of 5 seconds.

based on a time increment of 10 minutes.

5. Conclusions

D

Figure 6: Model concentration and temperature profiles.

320

TE

A model for the description of the drying kinetics of thin slabs is presented in this investigation. The resulting system of coupled partial differential equa-

EP

tions has been solved for two boundary conditions commonly encountered in the drying of foodstuff. The model has been validated by its application to three different experimental datasets. Aside from the good agreement found between the

AC C

model and the experimental data, the analysis has been focused on understand-

325

ing the model behaviour in relation with three parameters of influence, namely the effective diffusivity, convective coefficient and relative humidity. A convective and diffusivity controlled drying stages have been identified and explained by mean of parametrical variation studies of the aforementioned parameters. Furthermore, this method has been proposed as a substitute for the Biot num-

330

ber in the assessment whether a simplified exponential decay model is suitable for a particular drying experiment. 23

ACCEPTED MANUSCRIPT

Acknowledgement

RI PT

Thanks to Prof. Jos´e Mar´ıa Sa´ız Jabardo (Universidade de S˜ao Paulo) for reviewing the original manuscript. The funding provided by “Direcci´on de Inves335

tigaci´ on y Desarrollo-UTA” through grant number 2441-CU-P-2015 is gratefully acknowledge.

SC

Nomenclature Variables:

340

cp specific heat [J kg−1 K] C water concentration [kg m−3 ]

M AN U

A area [m2 ]

Def f effective diffusivity [m2 s−1 ]

k thermal conductivity [W m−1 K−1 ]

hd mass convection coefficient [m s−1 ] 345

hk heat convection coefficient [W m−2 K−1 ]

D

hl liquid enthalpy [J kg−1 ]

hlv liquid to vapor enthalpy [J kg−1 ]

TE

L sample thickness [m] Le Lewis number m mass [kg]

EP

350

N u Nusselt number P r Prandtl number

AC C

Re Reynolds number

RH relative humidity

355

x longitudinal coordinate [m] t time [s]

T temperature [K] Ai , Bi , Ci , Dij , Ji , Ki , Kij series expansion coefficient [K] Ej , Fji , Gji , Hjik , K, Oi series expansion coefficient [kg m−3 ]

360

Ix , i integral coefficient 24

ACCEPTED MANUSCRIPT

Mi series expansion coefficient

U, V, W, Z used in the method of separation of variables

365

Greek symbols: α thermal diffusivity [m s−2 ] ρ density [kg m−3 ]

SC

β model parameter [kgm−2 s−1 ] γ model parameter [kgm−2 s−1 K−1 ] Θ model parameter [K] Λ model parameter Υ model parameter [W m−1 ]

M AN U

370

ζ, υ, ξ, ϑ coefficient for series expansion [ m−1 ]

375

Subscripts: 0 initial condition

380

max maximum

EP

i, j, k cell index

TE

min minimum

D

eq equilibrium, related to water activity L evaluated at “x = L”

t time step

AC C

∞ distant field boundary condition 385

Superscripts:

dew dew point stage transition ∗ boundary at Ceq stage transition Acronyms:

390

RI PT

Q series expansion coefficient [kgm−5 ]

MABE mean absolute error RH relative humidity 25

ACCEPTED MANUSCRIPT

RI PT

References

ASHRAE, 2010. HandBook of Refrigeration. American Society of Heating, Re395

frigerating and Air-Conditioning Engineers, Atlanta, Georgia, United States.

Barati, E., Esfahani, J., 2013. A novel approach to evaluate the temperature

SC

during drying of food products with negligible external resistance to mass transfer. Journal of Food Engineering 114, 39–46.

400

M AN U

Bell, I.H., Wronski, J., Quoilin, S., Lemort, V., 2014.

Pure and

pseudo-pure fluid thermophysical property evaluation and the opensource

thermophysical

Engineering

Chemistry

property

library

Research

53,

coolprop.

2498–2508.

//pubs.acs.org/doi/abs/10.1021/ie4033999,

Industrial URL:

&

http:

doi:10.1021/ie4033999,

arXiv:http://pubs.acs.org/doi/pdf/10.1021/ie4033999. Belleghem, M.V., Steeman, M., Janssen, H., Janssens, A., Paepe, M.D., 2014.

D

405

Validation of a coupled heat, vapour and liquid moisture transport model for

TE

porous materials implemented in cfd. Building and Environment 81, 340–353. Bezerra, C.V., da Silva, L.H.M., Correa, D.F., Rodrigues, M., 2015. A modeling

410

EP

study for moisture diffusivities and moisture transfer coefficients in drying of passion fruit peel. International Journal of Heat and Mass Transfer 85, 750– 755.

AC C

Chiang, W., Petersen, J.N., 1987. Experimental measurement of temperature and moisture profiles during apple drying. Drying Technology 5, 25–49.

Choi, Y., Okos, M., 1986. Effects of temperature and composition on the thermal

415

properties of foods. 1 ed., M. LeMaguer and P. Jelen, eds. Elsevier Applied Science, London, United Kingdom.

Crank, J., 1975. The Mathematics of Diffusion. 2 ed., Clarendon Press, Oxford, United Kingdom. 26

ACCEPTED MANUSCRIPT

Datta, A., 2007. Porous media approaches to studying simultaneous heat and 420

mass transfer in food processes. i: Problem formulations. Journal of Food

RI PT

Engineering 80, 80–95.

Garc´ıa, J., 2018. C++ program for drying kinetics. Mendeley Data. http: http://dx.doi.org/10.17632/gfdkjcc8d9.1.

Garc´ıa-Alvarado, M., Pacheco-Aguirre, F.M., Ruiz-L´opez, I., 2014. Analytical solution of simultaneous heat and mass transfer equations during food drying. Journal of Food Engineering 142, 39–45.

SC

425

M AN U

GSL, 2017. GNU Scientific Library, v 1.16. https://www.gnu.org/software/ gsl/.

Hickey, C.J., Raspet, R., Slaton, W.V., 2000. Effects of thermal diffusion on 430

sound attenuation in evaporating and condensing gas-vapor mixtures in tubes. The Journal of the Acoustical Society of America 107, 1126–1130. Ilic, M., Turner, I., 1989. Convective drying of a consolidated slab of wet porous

D

material. International Journal of Heat and Mass Transfer 32, 2351–2362.

435

TE

Incropera, F.P., DeWitt, D.P., 1996. Fundamentals of heat and mass transfer. 4 ed., Wiley, John and Sons, New York, United States.

EP

Johnson, P.N.T., Brennan, J., 2000. Moisture sorption isotherm characteristics of plantain (musa aab). Journal of Food Engineering 44, 79–84.

AC C

Johnson, P.N.T., Brennan, J., Addo-Yobo, F., 1998. Air-drying characteristics of plantain (musa aab). Journal of Food Engineering 37, 233–242.

440

Kaya, A., Aydin, O., Demirtas, C., Akgn, M., 2007. An experimental study on the drying kinetics of quince. Desalination 212, 328–343.

Liu, J.Y., Cheng, S., 1991. Solutions of luikov equations of heat and mass transfer in capillary-porous bodies. International Journal of Heat and Mass Transfer 34, 1747–1754.

27

ACCEPTED MANUSCRIPT

445

Luikov, A.V., 1966. Heat and Mass Transfer in Capillary Porous Bodies. 1 ed.,

RI PT

Pergamon Press, New York. Ruiz-L´ opez, I., Ruiz-Espinosa, H., Arellanes-Lozada, P., B´arcenas-Pozos, M.,

Garc´ıa-Alvarado, M., 2012. Analytical model for variable moisture diffusivity estimation and drying simulation of shrinkable food products. Journal of Food 450

Engineering 108, 427–435.

SC

Tzempelikos, D.A., Mitrakos, D., Vouros, A.P., Bardakas, A.V., Filios, A.E., Margaris, D.P., 2015a. Numerical modelling of heat and mass transfer during

M AN U

convective drying of cylindrical quince slices. Journal of Food Engineering 156, 10–21. 455

Tzempelikos, D.A., Vouros, A.P., Bardakas, A.V., Filios, A.E., Margaris, D.P., 2015b. Experimental study on convective drying of quince slices and evaluation of thin-layer drying models. Engineering in Agriculture, Environment and Food 8, 169–177.

drying. Advances in heat transfer 31, 1–104.

TE

460

D

Whitaker, S., 1998. Coupled transport in multiphase systems: A theory of

White, F.M., 1995. Fluid Mechanics, 4th Ed. Mc Graw Hill, Essex, United Kingdom.

EP

Zhang, J., Datta, A.K., 2004. Some considerations in modeling of moisture transport in heating of hygroscopic materials. Drying Technology: An International Journal 22, 1983–2008.

AC C

465

Zhu, A., Shen, X., 2014. The model and mass transfer characteristics of convection drying of peach slices. International Journal of Heat and Mass Transfer 72, 345–351.

28

ACCEPTED MANUSCRIPT

RI PT

Annex A: Coefficients for the analytical solution

Table 8: Coefficients for case “a”.

2 Υ L XΥ · + · k 2 k L · ζi2 i=0

(35)

 2 (−1)i · e(−α·ζi ·t) −

1st stage

L + 6 · Def f

+ C0 + (γ + β · T (0)) ·

β·

∞ X i=0

L 6 · Def f

X Dik Ai Bi − − α · ζi2 α · ϑ2i α · ζk2

∞

Z

T (t) · dt =Θ · t +

Oi =

D

(52)

∞  X 2 Ai − · e(−αζi ·t) α · ζi2 i=0

AC C Q=−

hk · L · Λ k

1 L − · sin(2 · ϑi · L) 2 4 · ϑi 1 I2,i · L = − · (cos(ϑi · L) − 1) ϑi L 1 I3,i · L2 = − · cos(ϑi · L) + 2 · sin(ϑi · L) ϑi ϑi    L2 2 2 2·L I4,i · L3 = − + 3 · cos(ϑi · L) − 3 + 2 · sin(ϑi · L) ϑi ϑi ϑi ϑi I3,i − I4,i I5,i = I1,i I2,i I6,i = I1,i I1,i · L =

(40) (41)

(42) (43) (44) (45) (46) (47)

  Cmin − Cmax γ = hd · Cmax − C∞ − · T∞ T0 − T∞ Cmin − Cmax β = hd · T0 − T∞ j·π ξj = L Θ = T∞ −

Υ · k

(54) (55) (56)

! ∞ L X 2 + · 2 3 · sin(ζi · L) 2 L · ζi i=0

(57)

2 Υ · · (−1)i · sin(ζi · L) Ai = k L · ζi2

(58)

Bi = Ji · sin(ϑi · L)

(59)

Ci = Ki · sin(ϑi · L)

(60)

Dij = Kij · sin(ϑi · L)

(61)

(53)

Not computed

Not computed

γ + β · T (t∗ ) 2 · L · Def f

(62)

  h(t∗ ) − Ceq 2 (−1)i · (−1)i + Q · L2 − 2 · + ζi ζi ζi  !! ∞ X −ζi  + gj (t∗ ) · (−1)(i+j) · ξj2 − ζi2 j=1

2 · L

(50)

(51)

 ∞ X 2 2 Dij Bi · e(−αϑi ·t) + Ci · t − · e(−α·ζj ·t)  2 2 α · ϑi α · ζj j=0

T

−

− ϑi · L · tan−1 (ϑi · L) =

1 L

(49)

2 (−1)j · L ξj2 · Def f α · ζk2 − Def f · ξj2

β · α · ζk2 · Dik ·

·

!

TE

2 (−1)j · L ξj2 · Def f α · ϑ2i − Def f · ξj2

Hjik =

!!

k=0

2 (−1)j · L ξj2 · Def f α · ζi2 − Def f · ξj2

Gji =

(39)

k=0

X X 2 (−1)j Hjik Fji + Gji + − · L ξj2 · Def f i=0

β · α · ϑ2i · Bi ·

π · (2 · i + 1) 2·L

(48)

β · α · ζi2 · Ai ·

Fji =

(38)

EP

Concentration

(36)

(37)

∞

−

∞

Ej = (γ + β · T (0)) ·

Concentration

· sin(ζi · L)


2 (−1)j · sin(ζj · L) hk · Λ ·Υ· −I5,i · L2 · ζj2 + 2 · I6,i · 2 L2 · k 2 (ϑ2i − ζj2 ) ζj

K = − β · T (0) ·

Case “a”



∞ X I2,i − Ki − Kij Ji = (T0 − T∞ ) · I1,i j=0   ∞ 2 I2,i  L X 2 hk · Λ  · Υ · · · − − · sin(ζ · L) Ki = − j L · k2 ϑ2i I1,i 2 j=0 ζj3 · L2

Kij = −

2nd stage

1 ζi · L

ζi =

M AN U

Temperature

∞

VL = −

(34)

SC

hlv Cmin − Cmax · T0 − T∞ ρ · cp · Le(1−n) hlv · (C Υ = hk · max − C∞ ) ρ · cp · Le(1−n) Λ=1+

h(t) = β · T (t) ·



L 6 · Def f 2

  Z 1 − γ · t + β · T (t) · dt · + K L

gj (t) =Ej · e(−Def f ·ξj ·t) + (63)

+

∞ X k=0

29

∞  X i=0

2

(64)

2

Fji · e(−α·ζi ·t) + Gji · e(−α·ϑi ·t)

! 2

Hjik · e(−α·ζk ·t)

(65)

RI PT

ACCEPTED MANUSCRIPT

Table 9: Coefficients for case “c”.

hlv Cmin − Cmax · T0 − T∞ ρ · cp · Le(1−n) hlv · (C Υ = hk · max − C∞ ) ρ · cp · Le(1−n)

(67)

∞ X

Υ L Υ α Υ 1 2 (−α·ξj2 ·t) · − · ·t+ · · ·e k 3 k L k ξj2 L j=1

Ji = (T0 − T∞ ) ·

∞ X I2,i − Ki − Kij I1,i j=0

(69)

1 Υ 1 I2,i · · · ϑ2i k L I1,i Υ 2 I2,i 1 · · · Kij = − 2 ϑi − ξj2 k L I1,i Ki = −

K = − β · T (0) ·

+ C0 + (γ + β · T (0)) ·

1st stage

(70)

β·

∞ X i=0

∞

−

L 6 · Def f

X Dik Bi − α · ϑ2i α · ξk2 k=0

∞

j

2 (−1) β · α · · Dik · · 2 L ξj · Def f = 2 α · ξk − Def f · ξj2

!

(77)

(78)

(79)

1 L + · sin(2 · ϑi · L) 2 4 · ϑi 1 I2,i · L = · sin(ϑi · L) ϑi

I1,i · L =

(74) (75)

  Cmin − Cmax γ = hd · Cmax − C∞ − · T∞ T0 − T∞ Cmin − Cmax β = hd · T0 − T∞

(81) (82)

Bi = Ji · sin(ϑi · L)

(83)

Ci = Ki · sin(ϑi · L)

(84)

Dij = Kij · sin(ϑi · L)

(85) (86)

T (t) · dt = T∞ · t+   ∞ ∞ X X 2 2 Dij − Bi · e(−αϑi ·t) + Ci · t − · e(−α·ξj ·t)  2 2 α · ϑi α · ξj i=0 j=1

(80)

Not computed

T

(73)

1 L

EP

Z

·

TE

Hjik

!!

D

k=0

2 (−1)j · L ξj2 · Def f α · ϑ2i − Def f · ξj2

β · α · ϑ2i · Bi ·

ξk2

(72)

hk · L · Λ ϑi · L · tan(ϑi · L) = k

(76)

X X 2 (−1)j Hjik Ej = (γ + β · T (0)) · · 2 Gji + − L ξj · Def f i=0 Gji =

π·i L

(71)

∞

Concentration

Case “a”

L + 6 · Def f

(68)

ξi =

SC

VL (t) = −

(66)

M AN U

Temperature

Λ=1+

Not computed

∗

γ + β · T (t ) 2 · L · Def f

(87)

AC C Concentration

2nd stage

Q=−

Oi =

(88)

  h(t∗ ) − Ceq 2 (−1)i · (−1)i + Q · L2 − 2 · + ζi ζi ζi  !! ∞ X −ζi  + gj (t∗ ) · (−1)(i+j) · ξj2 − ζi2 j=1

2 · L

h(t) = β · T (t) ·

  Z 1 L − γ · t + β · T (t) · dt · + K 6 · Def f L



2

gj (t) = Ej · e(−Def f ·ξj ·t) + (89)

30

∞ X i=0

2

Gji · e(−α·ϑi ·t) +

∞ X k=1

(90)

! 2

Hjik · e(−α·ξk ·t)

(91)

Mj = 4 ·

1 − cos(υj · L) 2 · υj · L − sin(2 · υj · L)

− υj · L · tan−1 (υj · L) =

(92)

C

1st st.

T

Table 10: Coefficients for case “b”.

T

st.

I6,i =

∞ X

2

j=0

2

dew

Mj · e(−α·υj ·t

)

·

I7,ij I1,i

(95)

1 sin ((υj + ϑi ) · L) sin ((υj − ϑi ) · L) − 2 υj + ϑi υj − ϑi

I7,ij · L = −

(96)

Not computed

Same as case “a”, 2nd stage of drying.

D

T

C

Same as case “a”, 1st stage of drying.

C

st. 3

rd



(93)

M AN U



hk · L k

SC

(94)

Same as case “a”, 1st stage of drying but:

nd

Case “b”

C = C0

RI PT

ACCEPTED MANUSCRIPT

TE

Mj = 4 ·

sin(υj · L) 2 · υj · L + sin(2 · υj · L)

C

1

st

st.

T

Table 11: Coefficients for case “d”.

(99)

EP

C = C0

(97)

∞ X j=0

2

dew

Mj · e(−α·υj ·t

)

· I7,ij

T

AC C

2

nd

st.

I2,i =

I7,ij · L =

(100)

L

Z

0

cos(υj · x) · cos(ϑi · x) · dx =

1 sin ((υj + ϑi ) · L) sin ((υj − ϑi ) · L) + 2 υj + ϑi υj − ϑi

Same as case “c”, 1

(101)

stage of drying.

T

3

st



Not computed

C

st.

C



rd

Case “d”

Same as case “c”, 1st stage of drying but:

Same as case “c”, 2nd stage of drying.

31

υj · L · tan(υj · L) =

hk · L k

(98)

ACCEPTED MANUSCRIPT Highlights:

EP

TE D

M AN U

SC

RI PT

Modeling of one dimensional drying kinetics Analytical solution for the coupled heat and mass transfer differential equations Influence of the diffusivity, convection coefficient and humidity into drying curves Convective controlled and diffusion controlled drying stages

AC C

• • • •