- Email: [email protected]
Analytical and numerical study of drying of tomato in non-shrinkage and shrinkage model
Available online at www.sciencedirect.com
ScienceDirect Mathematics and Computers in Simulation 166 (2019) 253–265 www.elsevier.com/locate/matcom
Original articles
Analytical and numerical study of drying of tomato in non-shrinkage and shrinkage model Ebrahim Azhdaria ,∗, Aram Emamib a
Department of Mathematics, Salman Farsi University of Kazerun, Fars, Kazerun, Iran b Department of Mathematics, Fasa University, Fasa, Fars, Iran
Received 18 February 2018; received in revised form 26 May 2019; accepted 27 May 2019 Available online 3 June 2019
Abstract In this paper the transport of heat and moisture during drying of tomato is studied by considering its viscoelastic properties. The phenomenon is described by a set of three coupled partial differential equations that take into account heat transfer, diffusion of moisture and stress driven diffusion. The stability behavior of the model in non-shrinkage case is studied. In order to solve numerically the model we consider finite difference method, IMEX(implicit–explicit) method. Numerical simulations for comparison of shrinkage and non-shrinkage models, show an influence of parameters in qualitative agreement with the expected physical behavior. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Heat transfer; Drying; Non-Fickian; Viscoelastic behavior; Numerical simulation
1. Mathematical model formulation Tomato is one of the most important vegetables for the food industry. Especially dried tomato products are used as important ingredients for pizza, various vegetables, spicy dishes, and so on [1,2]. Several authors have studied heat and mass transfer models to describe the entrance of heat into the food and diffusion of moisture outwards towards the surface of food during drying [4,11,12] as well as drying of tomato [6,15]. However to the best of our knowledge the studying, analytically and numerically, of the drying of tomato from the viscoelastic point of view has not yet been addressed. Because of viscoelastic properties of tomato [9,14] the transfer of heat during drying of tomato does not obey the heat equation. In order to consider the viscoelastic properties of tomato we should add another term to the heat equation describing dissipation effects. In this paper we study a mathematical model to predict the influence of the mechanical properties of the tomato – viscoelasticity – in the transfer and diffusion rate. The viscoelastic behavior of the tomato is described by the Maxwell fluid model [10]. We consider a slice of tomato, as a homogeneous slice, in one dimension of thickness L, initially at a uniform temperature T0 and uniform moisture content M0 . The meaning and units of all variables are presented in Table 1. During air drying, heat is transferred from air towards the product center. Moisture diffuses outwards towards the ∗
Corresponding author. E-mail addresses: [email protected] (E. Azhdari), [email protected] (A. Emami).
https://doi.org/10.1016/j.matcom.2019.05.013 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.
254
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Table 1 Definition and units of the variables, parameters and constants. Symbol
Definition (unit)
Symbol
Definition (unit)
M M M0 T T T0 Tair Tsur σ σ σ0 D Dv µ∗ h Hv α µ k
Local moisture (w.b) Non-dimensional moisture Initial moisture (w.b) Local temperature (◦ C) Non-dimensional temperature Initial temperature (◦ C) Air temperature (◦ C) Surface temperature(◦ C) Stress (Pa) Non-dimensional stress Initial stress (Pa) Effective diffusion coefficient (m2 /s) Stress diffusion coefficient (m2 ◦ C/Pa s) Viscosity coefficient (Pa s) Heat transfer coefficient (W/m2 ◦ C) Heat of vaporization (kJ/kg) Thermal diffusivity (m2 /s) Viscosity of tomato (Pa s) Thermal conductivity (W/m ◦ C)
k1 k¯ E A t x L L0 Ω ∂Ω Ω Ωh ∂Ωh D2,h Dx D−x D−t Ih
Removal rate constant (w.b/s ◦ C) Constant (1/◦ C) Young modulus of the material (Pa) Constant (1/◦ C) Time (s) Space (m) Thickness of sample (m) Initial thickness (m) Domain of slice Boundary of domain Domain of slice and boundary Mesh nodes in domain Mesh nodes on the boundary Second order finite difference operator Forward finite difference (space) Backward finite difference (space) Backward finite difference (time) Grid points of domain
surface. This phenomenon is described by a system of partial differential equations (PDEs) as follows [11,12]: ∂M ∂2 M = D 2 − k1 T ; 0 < x < L , ∂t ∂x
(1)
∂ 2σ ∂T ∂2T = α 2 + Dv 2 ; ∂t ∂x ∂x
(2)
0 < x < L,
∂σ E ∂T + σ = −E ; 0 < x < L. ∂t µ(T ) ∂t Initial and boundary conditions are at t = 0: M = M0 , T = T0 , σ = σ0 , L = L 0 ,
(3)
(4)
at x = 0: ∂M ∂T = 0 and = 0, (5) ∂x ∂x at x = L: ∂T ∂σ ∂M M = 0 and − k − Dv + Hv D = h(Tsur − Tair ). (6) ∂x ∂x ∂x In Eqs. (1)–(6), M is local moisture content, T is local temperature, σ is stress, D represents the effective diffusion coefficient, Dv is stress diffusion coefficient, α is thermal diffusivity, E is the Young modulus of the material, µ is its viscosity and t is the time. k is thermal conductivity, Hv is the heat of vaporization, Tair is air temperature, Tsur is the surface temperature of the tomato and h is the heat transfer coefficient. The term −k1 T describes the decrease of the moisture inside the slice and the positive constant k1 represents the removal rate. 2 The viscoelastic behavior of the tomato is taken into account by the term Dv ∂∂ xσ2 (Eq. (2)). In fact σ is the stress response of the tomato to the strain exerted by the incoming heat. This term represents the opposition of the tomato to the entrance of the heat. Eq. (3) defines the viscoelastic behavior of the tomato as described by the Maxwell fluid model [10] ∂σ E ∂ϵ + σ =E , (7) ∂t µ ∂t
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
255
where ϵ is the strain produced by the incoming heat. As the tomato acts as a barrier to the entrance of the heat, then σ and ϵ are of opposite sign, and a minus sign should be considered in the right hand side of (7). We assume ¯ . Note that through this that the strain and the temperature are proportional, that is, there is k¯ > 0 such that ϵ = kT paper k¯ = 1. In this paper we take the Young modulus, E, as a constant but the viscosity, µ, depends on the temperature as follows [7]: µ = µ∗ e−AT , where µ∗ is the coefficient of kinematic viscosity and A is constant. The paper is organized as follows: In Section 2 the stability behavior of the model is presented. A finite difference method which solves the model numerically is described in Section 3. In Section 4 some numerical simulations are included with the comparison of shrinkage and non-shrinkage model and the study of the influence of some parameters. Finally conclusions are drawn in Section 5.
2. Stability analysis To simplify the presentation, we assume in this section that µ is constant. In order to study the stability behavior of Eqs. (1)–(3), since the initial temperature is zero, from (3) we easily obtain ∫ E E 2 t − µE (t−s) T (s)ds. (8) σ (t) = σ0 e− µ t − E T (t) + e µ 0 We deduce from (2) and (8), by assuming that the initial stress, σ0 , is constant, ∫ ∂T ∂2T Dv E 2 t − µE (t−s) ∂ 2 T = (α − Dv E) 2 + e (s)ds. ∂t ∂x µ ∂x2 0 Multiplying (1) and (9) by u(t) and v(t), respectively, and using Eqs. (5) and (6) we obtain ) ( ( ) ∂M ∂ M ∂u , u(t) = − D , − (k1 T, u(t)), ∂t ∂x ∂x
(9)
(10)
and (
) ( ) ∂ T ∂v ∂T , v(t) = − (α − Dv E) , ∂t ∂x ∂x ) ( 2 ∫ t E ∂v Dv E ∂T − (s), (t) ds e− µ (t−s) µ ∂x ∂x 0 ( ) + h T (L) − Tair v(L).
(11)
By taking u(t) = M(t) and v(t) = T (t) in (10) and (11), respectively, since Tair is bigger than or equal to the temperature on the boundary, we easily obtain 2 2 1 d M(t) = −D ∂ M (t) − (k1 T, M(t)), (12) 2 dt ∂x and 2 2 1 d T (t) ≤ −(α − Dv E) ∂ T (t) 2 dt ∂x ( ) 2 ∫ t E Dv E ∂T ∂T − e− µ (t−s) (s), (t) ds, µ ∂x ∂x 0 where ∥.∥ represents the usual norm in L 2 which is induced by the corresponding inner product (., .).
(13)
256
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Summing up Eqs. (12) and (13) and using the Cauchy–Schwarz inequality we have 2 2 2 2 1 d M(t) + 1 d T (t) + D ∂ M (t) + (α − Dv E) ∂ T (t) 2 dt 2 dt ∂x ∂x 2 2 1 ≤ k12 β 2 M(t) + 4β 2 T (t) ∫ 2 2 D 2 E 4 t − µE (t−s) ∂ T + ϵ 2 ∂ T (t) , + v2 2 e (s)ds 4ϵ µ ∂x ∂x 0 where ϵ, β ̸= 0, consequently we deduce 2 2 2 2 1 d M(t) + 1 d T (t) + D ∂ M (t) + (α − Dv E − ϵ 2 ) ∂ T (t) 2 dt 2 dt ∂x ∂x 2 2 1 ≤ k12 β 2 M(t) + 4β 2 T (t) ∫ ∫ t ∂ T 2 Dv2 E 4 t −2 µE (t−s) + 2 2 e ds ∂ x (s) ds, 4ϵ µ 0 0 and then ∫ t ∫ t ∂ M 2 ∂ T 2 ds + 2(α − Dv E − ϵ 2 ) (s) (s) ∂x ∂ x ds 0 0 2 2 ∫ t ∫ t 1 2 2 ≤ 2k1 β T (s) ds M(s) ds + 2β 2 0 0 2 ∫ ∫ Dv2 E 3 t s ∂T + (ν) dνds + Q(0), 4ϵ 2 µ 0 0 ∂ x 2 2 where Q(t) = M(t) + T (t) . If ϵ 2 is such that Q(t) + 2D
α − Dv E − ϵ 2 > 0,
(14)
we obtain E M T (t) ≤ C ∗
∫
t
E M T (s)ds + 0
1 Q(0), min{1, 2D, 2(α − Dv E − ϵ 2 )}
where ∫ t ∫ t ∂ M 2 ∂ T 2 E M T (t) = Q(t) + ∂ x (s) ds + ∂ x (s) ds, 0 0
(15)
and max{2k12 β 2 , 2β1 2 ,
∗
C =
Dv2 E 3 } 4ϵ 2 µ
min{1, 2D, 2(α − Dv E − ϵ 2 )}
.
(16)
Finally Gronwall’s Lemma [8] leads to E M T (t) ≤
1 ∗ Q(0)eC t . min{1, 2D, 2(α − Dv E − ϵ 2 )}
We have then proved the following result for the stability behavior of the model. Theorem 2.1. Let M and T be solutions of the variational problems (10) and (11), respectively, if (14) holds then E M T (t) ≤
1 ∗ Q(0)eC t , 2 min{1, 2D, 2(α − Dv E − ϵ )}
where E M T and C ∗ are defined by Eqs. (15) and (16), respectively.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
257
2 2 ∫t ∫t The inequality in Theorem 2.1 establishes that Q(t), 0 ∂∂Mx (s) ds and 0 ∂∂Tx (s) ds are bounded for bounded intervals of time. By following [3] and [13], this inequality can be improved by eliminating the exponential factor in its right hand side. 3. Numerical method In this section, we analyze the model, numerically, with and without structural change (in terms of shrinkage), which can affect the accuracy of the prediction model. To simplify the equations for numerical analysis, the following scaling group is used, Dt M T − T0 σ x τ = 2, X = , M= , T = , σ = . (17) L0 M0 Tair − T0 σ0 L0 3.1. Without shrinkage effect First we consider that there is no product shrinkage. That is, L is constant (L = L 0 ). Eqs. (1)–(3), after using (17), transform to ∂M ∂2 M = − k2 T , 0 < X < 1, ∂τ ∂ X2
(18)
∂T ∂ 2σ ∂2T , 0 < X < 1, = Le 2 + Dvnon ∂τ ∂X ∂ X2
(19)
∂σ ∂T + E non σ = −µnon , 0 < X < 1, ∂τ ∂τ
(20)
where Le =
E L 20 α Dv σ0 E(Tair − T0 ) , Dvnon = , E non = , µnon = , D D(Tair − T0 ) Dµ σ0
k1 L 20 (Tair − T0 ) , M0 D where Le, Dvnon , E non , µnon and k2 represent Lewis number [11], viscoelastic diffusion coefficient with respect to the diffusivity of the moisture, the Young modulus with respect to the diffusivity of the moisture and viscosity of the material, Young modulus with respect to the initial stress and removal rate constant with respect to the diffusivity of the moisture, respectively. Initial conditions become k2 =
T = 0, M = 1, σ = 1, at τ = 0,
(21)
and boundary conditions ∂M ∂T = 0 and = 0, on X = 0, ∂X ∂X M = 0 and −
∂T ∂σ ∂M −χ +λ = Bi(T sur − 1), on X = 1, ∂X ∂X ∂X
(22)
(23)
where Dv σ0 Hv D M0 h L0 , λ= , Bi = , k(Tair − T0 ) k(Tair − T0 ) k where χ, λ and Bi represent, respectively, stress diffusion coefficient with respect to the thermal conductivity of tomato, heat of vaporization and effective diffusion coefficient with respect to the thermal conductivity of tomato and Biot number for heat transfer [11]. χ=
258
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
We introduce a finite difference method to solve (18)–(23). Let Ω be the line (0, 1). We fix ∆x > 0 and we define in Ω the grid Ih = {xi , i = 0, 1, . . . , N , x0 = 0, x N = 1, xi − xi−1 = ∆x, i = 1, 2, . . . , N }. By Ωh and ∂Ωh we represent the mesh nodes of Ω h that are in Ω and on the boundary ∂Ω , respectively. In [0, t f ] we consider the following time grid {tn , n = 0, 1, . . . , N , t0 = 0, t N = t f , tn − tn−1 = ∆t, n = 1, 2, . . . , N }. Let phn (xi ) stands for an approximation of p(xi , tn ). To solve numerically the initial boundary value problems (18)–(23) we consider the IMEX (implicit–explicit) method defined by n+1
D−t M i
n+1
D−t T i
n+1
= D2,h M i
n
− k2 T i ; in Ω h
n+1
= LeD2,h T i
(24)
+ Dvnon D2,h σ in ; in Ω h n+1
D−t σ in+1 + E non σ in = −µnon D−t T i
(25)
; in Ω h
(26)
in Ωh : 0
0
M h = 1, T h = 0, σ 0h = 1, L 0 = 1,
(27)
on ∂Ωh (left): n+1
Dx M 0
n+1
= 0 and Dx T 0
= 0,
(28)
on ∂Ωh (right): n+1
MN
n+1
= 0 and − D−x T N
n+1
− χ D−x σ nN + λD−x M N
n+1
= Bi(T N
− 1),
(29)
for i = 1, . . . , N − 1 and n = 0, . . . , N − 1, where D−t , D−x and Dx denote, respectively, the backward finite difference operator with respect to the t-variable, the backward finite difference operator with respect to the x-variable and the forward finite difference operator with respect to the x-variable, D−t u n (xi ) =
u n (xi ) − u n−1 (xi ) , n = 1, . . . , N , ∆t
D−x u n (xi ) =
u n (xi ) − u n (xi−1 ) , i = 1, . . . , N , ∆x
u n (xi+1 ) − u n (xi ) , i = 1, . . . , N − 1, ∆x is the second order finite difference operator
Dx u n (xi ) = and D2,h
D2,h u(xi ) =
u(xi−1 ) − 2u(xi ) + u(xi+1 ) , i = 1, . . . , N − 1. (∆x)2
3.2. With shrinkage effect The shape and size of the sample particles are constantly changing as a result of water removal during drying. That is L is changing with time (L = L(t)). Using the same non-dimensional scale variable, (17), and new variable l(τ ) = L(t) , the domain of Eqs. (18)–(23) change to the interval 0 < X < l(τ ). By using the L0 transformationξ = l(τX ) [5], the surface interface corresponds to a fixed value ξ = 1 and Eqs. (18)–(23) become
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
259
1 ∂2 M ξ dl ∂ M ∂M = 2 + − k2 T , 0 < ξ < 1, 2 ∂τ l (τ ) ∂ξ l(τ ) dτ ∂ξ
(30)
Le ∂ 2 T Dvnon ∂ 2 σ ξ dl ∂ T ∂T = 2 + + , 0 < ξ < 1, 2 2 2 ∂τ l (τ ) ∂ξ l (τ ) ∂ξ l(τ ) dτ ∂ξ
(31)
( ) ξ dl ∂T ∂σ ∂T ∂σ µnon − µnon + E non σ = + , 0 < ξ < 1. ∂τ l(τ ) dτ ∂ξ ∂ξ ∂τ
(32)
M(ξ, 0) = 1, T (ξ, 0) = 0, σ (ξ, 0) = 1,
(33)
∂T ∂M = 0 and = 0, ξ = 0, ∂ξ ∂ξ
(34)
∂σ ∂M ∂T −χ +λ = l(τ )Bi(T sur − 1), ξ = 1. (35) ∂ξ ∂ξ ∂ξ The further condition is needed on the moving interface to determine the position of the interface itself (see [5]). In this model it is M = 0 and −
dl M0 ∂ M = , ξ = 1, (36) dτ l(τ ) ∂ξ with l(0) = 1. Eq. (36) is known as the ‘Stefan condition’ and it expresses the heat balance on the interface. Eqs. (30)–(32) were solved together with (36), which involves a shrinkage equation obtained using the IMEX method as above. Discretizing Eq. (30)–(32) with a central difference approximation gives, ( n+1 n+1 n+1 ) n+1 n M i−1 − 2M i + M i+1 Mi − Mi 1 = 2 ∆τ l (n + 1) (∆ξ )2 ( )( n+1 n+1 ) M i+1 − M i−1 ξi l(n + 1) − l(n) n + (37) − k2 T i , l(n + 1) ∆τ 2∆ξ ( n+1 n+1 n+1 ) n T i−1 − 2T i + T i+1 − Ti Le = 2 ∆τ l (n + 1) (∆ξ )2 ( n n ) σ i−1 − 2σ in + σ i+1 Dvnon + 2 l (n + 1) (∆ξ )2 ( )( n+1 n+1 ) T i+1 − T i−1 l(n + 1) − l(n) ξi , + l(n + 1) ∆τ 2∆ξ
n+1
Ti
(38)
( )( n+1 n+1 n ) σ n − σ i−1 T − T i−1 σ in+1 − σ in ξi l(n + 1) − l(n) = µnon i+1 + i+1 ∆τ l(n + 1) ∆τ 2∆ξ 2∆ξ n+1
−
E non σin
− µnon
Ti
n
− Ti , ∆τ
(39)
respectively, for i = 1, . . . , N − 1 and n = 0, . . . , N − 1. Discretization of Eq. (36) by using the method proposed by Furzeland in [5] we have ( ) n n n 3M N − 4M N −1 + M N −2 l(n + 1) − l(n) M0 = , (40) ∆τ l(n + 1) 2∆ξ for n = 0, . . . , N − 1.
260
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265 Table 2 Input parameters used in the simulations of drying. Parameter
Value
D α Dv E µ∗ Hv Tair k1 k
8 × 10−10
m2 /s
4 × 10−9 m2 /s 8 × 10−12 (m2 ◦ C)/(Pa s) 1 × 10−9 Pa 1 × 10−4 Pa s 2345 kJ/kg 60 ◦ C 1 × 10−10 (w.b)/(s ◦ C) 0.475 W/m ◦ C
Parameter
Value
h M0 T0 σ0 L0 ∆x ∆t A k¯
30 W/m2 ◦ C 0.8 w.b 0◦ C 5 × 10−2 Pa 5 × 10−3 m 2 × 10−2 2 × 10−7 1 × 10−4 (1/◦ C) 1(1/◦ C)
Table 3 Errors and convergence orders for the temperature and moisture for different step-size ∆x. ∆x
E∆x (T )
R(T )
E∆x (M)
R(M)
0.01 0.005 0.004 0.002
3.2 × 10−3 1.6 × 10−3 1.3 × 10−3 5.1303 × 10−4
– 1 0.9305 1.3414
1.57 × 10−2 1.76 × 10−2 1.64 × 10−2 9.7 × 10−3
– −0.1648 0.3165 0.7576
Table 4 Errors and convergence orders for the temperature and moisture for different step-size ∆t. ∆t
E∆t (T )
2 × 10−3
3.68 × 10−2
2 × 10−4 2 × 10−5 2 × 10−6
8.2 × 10−3 1 × 10−3 9.9433 × 10−5
R(T )
E∆t (M)
R(M)
– 0.6520 0.9138 1.0025
2.9629 × 10−6
– 0.7874 0.9143 1.0274
4.8337 × 10−7 5.8877 × 10−8 5.5276 × 10−9
4. Numerical simulations In this section we obtain some numerical simulations for Eqs. (1)–(6) without shrinkage effect and with shrinkage effect by using Eqs. (24)–(29) and (37)–(40), respectively. For drying simulation, the input parameters are given in Table 2 for generic drying conditions for fruit [12] and parameters which are specially used for viscoelastic behavior of tomato. We start by analyzing the numerical convergence properties of the method (37)–(40). The behavior of the error E∆x and convergence order R are illustrated for different step-size ∆x when the time step-size is fixed at time t f = 1 × 10−3 . The error was computed using a reference solution obtained with ∆x = 0.001 and ∆t = 2 × 10−7 for T and M. In Table 3, E∆x is defined by ⎛ ⎞1 2 ∑ n n n n 2⎠ ⎝ E∆x (c) = ∥c∆x − c¯∆x ∥ = (c∆x ( p) − c¯∆x ( p)) , (41) p∈Ω
where
n c¯∆x
is the reference solution, and the rate R is defined by
R(c) =
E∆x1 (c) ) ∆x2 (c)
ln( E
1 ln( ∆x ) ∆x
,
(42)
2
where ∆x1 and ∆x2 represent two consecutive step sizes. We note that the numerical results presented in Table 3, suggest the convergence of the method (37)–(40) with first order convergence rate in space. Table 4 illustrates the behavior of the error E∆t and convergence order R, replacing ∆x by ∆t in Eqs. (41)– (42), for different time step-sizes when the step-size ∆x is fixed at time t f = 0.5. The error was computed using a reference solution obtained with ∆x = 0.02 and ∆t = 2 × 10−7 for T and M.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
261
Fig. 1. Profile of moisture through the sample with elapsed time τ = 0 − 1 in step ∆t = 2 × 10−7 for (a) Non-shrinkage model (b) shrinkage model (c) location of the surface.
Fig. 1(a, b) shows the moisture profile through the sample of the tomato with increasing time for the shrinkage and non-shrinkage models. Moisture decreased as time increased, but this process was a little slower for the nonshrinkage model compared to the shrinkage model. Fig. 1(c) shows the location of the free surface during shrinkage; the shrinkage of the product is shown by the change in the thickness of the tomato. The final value of shrinkage in this model can be determined readily as only solid mass is left at the end of drying, which in this case is 0.4576. Fig. 2 shows a temperature profile in the sample as a function of both position and time during drying for the non-shrinkage and the shrinkage model. The temperature at each location increases with the drying time, which is due to the higher drying temperature. Due to the difference between the air temperature and the sample temperature, the temperature profile rises rapidly in the early period of heating. As the heating period progresses, the rise in temperature attains an almost uniform profile. The overall dimension varies with time depending on the amount of moisture loss for each section, which can be seen in Fig. 2(b). The temperature distribution increased much faster for the shrinkage model compared with the non-shrinkage model. Fig. 3 represents the comparison of moisture and temperature gradients at the center of the sample, with and without the shrinkage effect. It can be seen that the temperature with shrinkage predicts a higher temperature at a specific time, compared to that with non-shrinkage. However, the predicted moisture content decreased faster when the shrinkage effect was included in the model. As a result of the thickness of the sample and decreases due to shrinkage, the moisture has less distance to cover and hence it reaches the surface faster, with the result that the time taken to dry a sample becomes shorter. Therefore, the model with shrinkage needs a shorter time for drying. Fig. 4 shows the profile of moisture in the sample at time τ = 0.25 in step of 2 × 10−7 , when different air temperatures are used during drying. An increase in drying temperature can influence the drying kinetics that can cause a decrease in the moisture of sample. We plot, in Fig. 5, the dependence of temperature at the center on the Dvnon which is proportional to the viscoelastic diffusion coefficient, Dv . In this figure we observe that the sample is acting as a barrier to the entrance
262
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Fig. 2. Profile of temperature through the sample with elapsed time τ = 0 − 1 in step 2 × 10−7 for (a) Non-shrinkage model (b) shrinkage model.
Fig. 3. Comparison of moisture and temperature ratio at the center for shrinkage and non-shrinkage model.
Fig. 4. Comparison of moisture at time τ = 0.25 in step of 2 × 10−7 with different air temperatures (for shrinkage model).
of heat. This is because of the viscoelastic properties of the tomato. According to this description, an increase in Dv the temperature at the center of the sample is decreasing and it causes the moisture not decreasing. The detection of this region is a very important aspect in terms of food safety. Fig. 6 shows that by increasing the value of Le, the temperature of the center increases. This is consistent with the result of [12], that is, for the value of Le more than 100, the interior temperature is approximately equal to the air temperature. Fig. 7 shows that by increasing the value of Bi, the temperature increases rapidly.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
263
Fig. 5. Influence of the Dv , viscoelastic diffusion coefficient, on the temperature at the center of the sample.
Fig. 6. Temperature at the center for different values of Le.
Fig. 7. Temperature at the surface for different values of Bi.
Fig. 8 shows the experiment measurement of shrinkage L/L 0 taken from experiment drying of tomato [6] at 60 ◦ C, which was compared with the simulation model of shrinkage. The thickness is assumed to decrease linearly
264
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Fig. 8. Comparison between experiment data with simulation study of variation of shrinkage with respect to the moisture content.
with a reduction in moisture content. The trend deviates with higher moisture content but there is good agreement at lower moisture content. 5. Conclusion A model to simulate the drying of tomato has been studied. The whole process is described by a set of partial differential equations that take into account the transfer of heat from air towards the product center, diffusion of moisture outwards towards the surface and stress driven diffusion. The stability of the model subject to nonshrinkage was established under a mathematical relation between the parameters involved in the equations. To solve the system of equations an IMEX method was proposed and its convergence was numerically studied. The numerical simulations show qualitative agreement with the physical expected behavior. Moisture decreases more quickly in the shrinkage model, which concludes that the model with shrinkage effect needs a shorter time for drying. The viscoelastic properties of the tomato are shown to be an effective control mechanism to delay the heat transfer. In future work we are planning to complete the model by considering three dimensionality and studying the stability analysis in time-varying domain. From an analytical point of view the convergence analysis of the IMEX methods used in the simulations will deserve our attention in the future. References [1] C.T. Akanbi, R.S. Adeyemi, A. Ojo, Drying characteristics and sorption isotherm of tomato slices, J. Food Eng. 73 (2006) 157–163, http://dx.doi.org/10.1016/j.jfoodeng.2005.01.015. [2] D. Arslan, M.M. Ozcan, Drying of tomato slices: changes in drying kinetics, mineral contents, antioxidant activity and color parameters, CytA-J. Food 9 (2011) 229–236, http://dx.doi.org/10.1080/19476337.2010.522734. [3] E. Azhdari, J.A. Ferreira, P. de Oliveira, P.M. da Silva, Diffusion, viscoelasticity and erosion: Analytical study and medical applications, J. Comput. Appl. Math. 275 (2015) 489–501, http://dx.doi.org/10.1016/j.cam.2014.01.025. [4] B. Broyart, G. Trystram, Modelling of heat and mass transfer phenomena and quality changes during continuous biscuit baking using both deductive and inductive (neural network modelling principles), Food Bioprod. Process. 81 (2003) 316–326, http://dx.doi.org/10. 1205/096030803322756402. [5] J. Crank, Free and Moving Boundary Problem, New York Clarendon-Oxford Science Publication, 1984. [6] B. Dianda, M. Ousmane, T.K. S. Kam, D.J. Bathiébo, Experimental study of the kinetics and shrinkage of tomato slices in convective drying, Afr. J. Food Sci. 9(5) (2015) 262–271, http://dx.doi.org/10.5897/AJFS2015.1298. [7] S.K. Ghosh, G.C. Shit, J.C. Misra, Heat transfer in hydromagnetic fluid flow: Study of temperature dependence of fluid viscosity, J. Appl. Fluid Mech. 7 (2014) 633–640. [8] T. Gronwall, Note on the derivative with respect to a parameter of the solutions of a system of differential equations, Ann. of Math. 20 (1919) 292–296, http://dx.doi.org/10.2307/1967124. [9] R. Jakman, A.G. Marangoni, D.W. Stanley, The effects of turgor pressure on puncture and viscoelastic properties of tomato tissue, J. Text. Stud. 23(4): (1992) 491–505, http://dx.doi.org/10.1111/j.1745-4603.1992.tb00036.x. [10] R.G. Owens, T.N. Phillips, Computational Rheology, Imperial College Press, 2002.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
265
[11] N.A. Shahari, Mathematical Modelling of Drying Food Products: Application to Tropical Fruits (Ph.D. Thesis), University of Nottingham, 2012. [12] N.A. Shahari, S. Hibberd, Analysis of single phase moisture and heat model of food drying, in: 4th International conference of Mathematical Models in Engineering and Computer science, 2013, pp. 137–142. [13] V. Thomée, B. Wahlbin, Long-time numerical solution of a parabolic equation with memory, Math. Comp. 62 (1994) 477–496. [14] I. Verlent, M. Hendickx, P. Rovere, P. Moldenaers, A.V. Loey, Rheological properties of tomato based products after thermal and high-pressure treatment, J. Food Sci. 71(3) (2006) 5243–5248, http://dx.doi.org/10.1111/j.1365-2621.2006.tb15648.x. [15] G. Xanthopoulos, S. Yanniotis, E. Talaiporou, Influence of salting on drying kinetics and water diffusivity of tomato halves, Int. J. Food Prop. 15(4) (2012) 847–863, http://dx.doi.org/10.1080/10942912.2010.506018.
ScienceDirect Mathematics and Computers in Simulation 166 (2019) 253–265 www.elsevier.com/locate/matcom
Original articles
Analytical and numerical study of drying of tomato in non-shrinkage and shrinkage model Ebrahim Azhdaria ,∗, Aram Emamib a
Department of Mathematics, Salman Farsi University of Kazerun, Fars, Kazerun, Iran b Department of Mathematics, Fasa University, Fasa, Fars, Iran
Received 18 February 2018; received in revised form 26 May 2019; accepted 27 May 2019 Available online 3 June 2019
Abstract In this paper the transport of heat and moisture during drying of tomato is studied by considering its viscoelastic properties. The phenomenon is described by a set of three coupled partial differential equations that take into account heat transfer, diffusion of moisture and stress driven diffusion. The stability behavior of the model in non-shrinkage case is studied. In order to solve numerically the model we consider finite difference method, IMEX(implicit–explicit) method. Numerical simulations for comparison of shrinkage and non-shrinkage models, show an influence of parameters in qualitative agreement with the expected physical behavior. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Heat transfer; Drying; Non-Fickian; Viscoelastic behavior; Numerical simulation
1. Mathematical model formulation Tomato is one of the most important vegetables for the food industry. Especially dried tomato products are used as important ingredients for pizza, various vegetables, spicy dishes, and so on [1,2]. Several authors have studied heat and mass transfer models to describe the entrance of heat into the food and diffusion of moisture outwards towards the surface of food during drying [4,11,12] as well as drying of tomato [6,15]. However to the best of our knowledge the studying, analytically and numerically, of the drying of tomato from the viscoelastic point of view has not yet been addressed. Because of viscoelastic properties of tomato [9,14] the transfer of heat during drying of tomato does not obey the heat equation. In order to consider the viscoelastic properties of tomato we should add another term to the heat equation describing dissipation effects. In this paper we study a mathematical model to predict the influence of the mechanical properties of the tomato – viscoelasticity – in the transfer and diffusion rate. The viscoelastic behavior of the tomato is described by the Maxwell fluid model [10]. We consider a slice of tomato, as a homogeneous slice, in one dimension of thickness L, initially at a uniform temperature T0 and uniform moisture content M0 . The meaning and units of all variables are presented in Table 1. During air drying, heat is transferred from air towards the product center. Moisture diffuses outwards towards the ∗
Corresponding author. E-mail addresses: [email protected] (E. Azhdari), [email protected] (A. Emami).
https://doi.org/10.1016/j.matcom.2019.05.013 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.
254
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Table 1 Definition and units of the variables, parameters and constants. Symbol
Definition (unit)
Symbol
Definition (unit)
M M M0 T T T0 Tair Tsur σ σ σ0 D Dv µ∗ h Hv α µ k
Local moisture (w.b) Non-dimensional moisture Initial moisture (w.b) Local temperature (◦ C) Non-dimensional temperature Initial temperature (◦ C) Air temperature (◦ C) Surface temperature(◦ C) Stress (Pa) Non-dimensional stress Initial stress (Pa) Effective diffusion coefficient (m2 /s) Stress diffusion coefficient (m2 ◦ C/Pa s) Viscosity coefficient (Pa s) Heat transfer coefficient (W/m2 ◦ C) Heat of vaporization (kJ/kg) Thermal diffusivity (m2 /s) Viscosity of tomato (Pa s) Thermal conductivity (W/m ◦ C)
k1 k¯ E A t x L L0 Ω ∂Ω Ω Ωh ∂Ωh D2,h Dx D−x D−t Ih
Removal rate constant (w.b/s ◦ C) Constant (1/◦ C) Young modulus of the material (Pa) Constant (1/◦ C) Time (s) Space (m) Thickness of sample (m) Initial thickness (m) Domain of slice Boundary of domain Domain of slice and boundary Mesh nodes in domain Mesh nodes on the boundary Second order finite difference operator Forward finite difference (space) Backward finite difference (space) Backward finite difference (time) Grid points of domain
surface. This phenomenon is described by a system of partial differential equations (PDEs) as follows [11,12]: ∂M ∂2 M = D 2 − k1 T ; 0 < x < L , ∂t ∂x
(1)
∂ 2σ ∂T ∂2T = α 2 + Dv 2 ; ∂t ∂x ∂x
(2)
0 < x < L,
∂σ E ∂T + σ = −E ; 0 < x < L. ∂t µ(T ) ∂t Initial and boundary conditions are at t = 0: M = M0 , T = T0 , σ = σ0 , L = L 0 ,
(3)
(4)
at x = 0: ∂M ∂T = 0 and = 0, (5) ∂x ∂x at x = L: ∂T ∂σ ∂M M = 0 and − k − Dv + Hv D = h(Tsur − Tair ). (6) ∂x ∂x ∂x In Eqs. (1)–(6), M is local moisture content, T is local temperature, σ is stress, D represents the effective diffusion coefficient, Dv is stress diffusion coefficient, α is thermal diffusivity, E is the Young modulus of the material, µ is its viscosity and t is the time. k is thermal conductivity, Hv is the heat of vaporization, Tair is air temperature, Tsur is the surface temperature of the tomato and h is the heat transfer coefficient. The term −k1 T describes the decrease of the moisture inside the slice and the positive constant k1 represents the removal rate. 2 The viscoelastic behavior of the tomato is taken into account by the term Dv ∂∂ xσ2 (Eq. (2)). In fact σ is the stress response of the tomato to the strain exerted by the incoming heat. This term represents the opposition of the tomato to the entrance of the heat. Eq. (3) defines the viscoelastic behavior of the tomato as described by the Maxwell fluid model [10] ∂σ E ∂ϵ + σ =E , (7) ∂t µ ∂t
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
255
where ϵ is the strain produced by the incoming heat. As the tomato acts as a barrier to the entrance of the heat, then σ and ϵ are of opposite sign, and a minus sign should be considered in the right hand side of (7). We assume ¯ . Note that through this that the strain and the temperature are proportional, that is, there is k¯ > 0 such that ϵ = kT paper k¯ = 1. In this paper we take the Young modulus, E, as a constant but the viscosity, µ, depends on the temperature as follows [7]: µ = µ∗ e−AT , where µ∗ is the coefficient of kinematic viscosity and A is constant. The paper is organized as follows: In Section 2 the stability behavior of the model is presented. A finite difference method which solves the model numerically is described in Section 3. In Section 4 some numerical simulations are included with the comparison of shrinkage and non-shrinkage model and the study of the influence of some parameters. Finally conclusions are drawn in Section 5.
2. Stability analysis To simplify the presentation, we assume in this section that µ is constant. In order to study the stability behavior of Eqs. (1)–(3), since the initial temperature is zero, from (3) we easily obtain ∫ E E 2 t − µE (t−s) T (s)ds. (8) σ (t) = σ0 e− µ t − E T (t) + e µ 0 We deduce from (2) and (8), by assuming that the initial stress, σ0 , is constant, ∫ ∂T ∂2T Dv E 2 t − µE (t−s) ∂ 2 T = (α − Dv E) 2 + e (s)ds. ∂t ∂x µ ∂x2 0 Multiplying (1) and (9) by u(t) and v(t), respectively, and using Eqs. (5) and (6) we obtain ) ( ( ) ∂M ∂ M ∂u , u(t) = − D , − (k1 T, u(t)), ∂t ∂x ∂x
(9)
(10)
and (
) ( ) ∂ T ∂v ∂T , v(t) = − (α − Dv E) , ∂t ∂x ∂x ) ( 2 ∫ t E ∂v Dv E ∂T − (s), (t) ds e− µ (t−s) µ ∂x ∂x 0 ( ) + h T (L) − Tair v(L).
(11)
By taking u(t) = M(t) and v(t) = T (t) in (10) and (11), respectively, since Tair is bigger than or equal to the temperature on the boundary, we easily obtain 2 2 1 d M(t) = −D ∂ M (t) − (k1 T, M(t)), (12) 2 dt ∂x and 2 2 1 d T (t) ≤ −(α − Dv E) ∂ T (t) 2 dt ∂x ( ) 2 ∫ t E Dv E ∂T ∂T − e− µ (t−s) (s), (t) ds, µ ∂x ∂x 0 where ∥.∥ represents the usual norm in L 2 which is induced by the corresponding inner product (., .).
(13)
256
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Summing up Eqs. (12) and (13) and using the Cauchy–Schwarz inequality we have 2 2 2 2 1 d M(t) + 1 d T (t) + D ∂ M (t) + (α − Dv E) ∂ T (t) 2 dt 2 dt ∂x ∂x 2 2 1 ≤ k12 β 2 M(t) + 4β 2 T (t) ∫ 2 2 D 2 E 4 t − µE (t−s) ∂ T + ϵ 2 ∂ T (t) , + v2 2 e (s)ds 4ϵ µ ∂x ∂x 0 where ϵ, β ̸= 0, consequently we deduce 2 2 2 2 1 d M(t) + 1 d T (t) + D ∂ M (t) + (α − Dv E − ϵ 2 ) ∂ T (t) 2 dt 2 dt ∂x ∂x 2 2 1 ≤ k12 β 2 M(t) + 4β 2 T (t) ∫ ∫ t ∂ T 2 Dv2 E 4 t −2 µE (t−s) + 2 2 e ds ∂ x (s) ds, 4ϵ µ 0 0 and then ∫ t ∫ t ∂ M 2 ∂ T 2 ds + 2(α − Dv E − ϵ 2 ) (s) (s) ∂x ∂ x ds 0 0 2 2 ∫ t ∫ t 1 2 2 ≤ 2k1 β T (s) ds M(s) ds + 2β 2 0 0 2 ∫ ∫ Dv2 E 3 t s ∂T + (ν) dνds + Q(0), 4ϵ 2 µ 0 0 ∂ x 2 2 where Q(t) = M(t) + T (t) . If ϵ 2 is such that Q(t) + 2D
α − Dv E − ϵ 2 > 0,
(14)
we obtain E M T (t) ≤ C ∗
∫
t
E M T (s)ds + 0
1 Q(0), min{1, 2D, 2(α − Dv E − ϵ 2 )}
where ∫ t ∫ t ∂ M 2 ∂ T 2 E M T (t) = Q(t) + ∂ x (s) ds + ∂ x (s) ds, 0 0
(15)
and max{2k12 β 2 , 2β1 2 ,
∗
C =
Dv2 E 3 } 4ϵ 2 µ
min{1, 2D, 2(α − Dv E − ϵ 2 )}
.
(16)
Finally Gronwall’s Lemma [8] leads to E M T (t) ≤
1 ∗ Q(0)eC t . min{1, 2D, 2(α − Dv E − ϵ 2 )}
We have then proved the following result for the stability behavior of the model. Theorem 2.1. Let M and T be solutions of the variational problems (10) and (11), respectively, if (14) holds then E M T (t) ≤
1 ∗ Q(0)eC t , 2 min{1, 2D, 2(α − Dv E − ϵ )}
where E M T and C ∗ are defined by Eqs. (15) and (16), respectively.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
257
2 2 ∫t ∫t The inequality in Theorem 2.1 establishes that Q(t), 0 ∂∂Mx (s) ds and 0 ∂∂Tx (s) ds are bounded for bounded intervals of time. By following [3] and [13], this inequality can be improved by eliminating the exponential factor in its right hand side. 3. Numerical method In this section, we analyze the model, numerically, with and without structural change (in terms of shrinkage), which can affect the accuracy of the prediction model. To simplify the equations for numerical analysis, the following scaling group is used, Dt M T − T0 σ x τ = 2, X = , M= , T = , σ = . (17) L0 M0 Tair − T0 σ0 L0 3.1. Without shrinkage effect First we consider that there is no product shrinkage. That is, L is constant (L = L 0 ). Eqs. (1)–(3), after using (17), transform to ∂M ∂2 M = − k2 T , 0 < X < 1, ∂τ ∂ X2
(18)
∂T ∂ 2σ ∂2T , 0 < X < 1, = Le 2 + Dvnon ∂τ ∂X ∂ X2
(19)
∂σ ∂T + E non σ = −µnon , 0 < X < 1, ∂τ ∂τ
(20)
where Le =
E L 20 α Dv σ0 E(Tair − T0 ) , Dvnon = , E non = , µnon = , D D(Tair − T0 ) Dµ σ0
k1 L 20 (Tair − T0 ) , M0 D where Le, Dvnon , E non , µnon and k2 represent Lewis number [11], viscoelastic diffusion coefficient with respect to the diffusivity of the moisture, the Young modulus with respect to the diffusivity of the moisture and viscosity of the material, Young modulus with respect to the initial stress and removal rate constant with respect to the diffusivity of the moisture, respectively. Initial conditions become k2 =
T = 0, M = 1, σ = 1, at τ = 0,
(21)
and boundary conditions ∂M ∂T = 0 and = 0, on X = 0, ∂X ∂X M = 0 and −
∂T ∂σ ∂M −χ +λ = Bi(T sur − 1), on X = 1, ∂X ∂X ∂X
(22)
(23)
where Dv σ0 Hv D M0 h L0 , λ= , Bi = , k(Tair − T0 ) k(Tair − T0 ) k where χ, λ and Bi represent, respectively, stress diffusion coefficient with respect to the thermal conductivity of tomato, heat of vaporization and effective diffusion coefficient with respect to the thermal conductivity of tomato and Biot number for heat transfer [11]. χ=
258
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
We introduce a finite difference method to solve (18)–(23). Let Ω be the line (0, 1). We fix ∆x > 0 and we define in Ω the grid Ih = {xi , i = 0, 1, . . . , N , x0 = 0, x N = 1, xi − xi−1 = ∆x, i = 1, 2, . . . , N }. By Ωh and ∂Ωh we represent the mesh nodes of Ω h that are in Ω and on the boundary ∂Ω , respectively. In [0, t f ] we consider the following time grid {tn , n = 0, 1, . . . , N , t0 = 0, t N = t f , tn − tn−1 = ∆t, n = 1, 2, . . . , N }. Let phn (xi ) stands for an approximation of p(xi , tn ). To solve numerically the initial boundary value problems (18)–(23) we consider the IMEX (implicit–explicit) method defined by n+1
D−t M i
n+1
D−t T i
n+1
= D2,h M i
n
− k2 T i ; in Ω h
n+1
= LeD2,h T i
(24)
+ Dvnon D2,h σ in ; in Ω h n+1
D−t σ in+1 + E non σ in = −µnon D−t T i
(25)
; in Ω h
(26)
in Ωh : 0
0
M h = 1, T h = 0, σ 0h = 1, L 0 = 1,
(27)
on ∂Ωh (left): n+1
Dx M 0
n+1
= 0 and Dx T 0
= 0,
(28)
on ∂Ωh (right): n+1
MN
n+1
= 0 and − D−x T N
n+1
− χ D−x σ nN + λD−x M N
n+1
= Bi(T N
− 1),
(29)
for i = 1, . . . , N − 1 and n = 0, . . . , N − 1, where D−t , D−x and Dx denote, respectively, the backward finite difference operator with respect to the t-variable, the backward finite difference operator with respect to the x-variable and the forward finite difference operator with respect to the x-variable, D−t u n (xi ) =
u n (xi ) − u n−1 (xi ) , n = 1, . . . , N , ∆t
D−x u n (xi ) =
u n (xi ) − u n (xi−1 ) , i = 1, . . . , N , ∆x
u n (xi+1 ) − u n (xi ) , i = 1, . . . , N − 1, ∆x is the second order finite difference operator
Dx u n (xi ) = and D2,h
D2,h u(xi ) =
u(xi−1 ) − 2u(xi ) + u(xi+1 ) , i = 1, . . . , N − 1. (∆x)2
3.2. With shrinkage effect The shape and size of the sample particles are constantly changing as a result of water removal during drying. That is L is changing with time (L = L(t)). Using the same non-dimensional scale variable, (17), and new variable l(τ ) = L(t) , the domain of Eqs. (18)–(23) change to the interval 0 < X < l(τ ). By using the L0 transformationξ = l(τX ) [5], the surface interface corresponds to a fixed value ξ = 1 and Eqs. (18)–(23) become
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
259
1 ∂2 M ξ dl ∂ M ∂M = 2 + − k2 T , 0 < ξ < 1, 2 ∂τ l (τ ) ∂ξ l(τ ) dτ ∂ξ
(30)
Le ∂ 2 T Dvnon ∂ 2 σ ξ dl ∂ T ∂T = 2 + + , 0 < ξ < 1, 2 2 2 ∂τ l (τ ) ∂ξ l (τ ) ∂ξ l(τ ) dτ ∂ξ
(31)
( ) ξ dl ∂T ∂σ ∂T ∂σ µnon − µnon + E non σ = + , 0 < ξ < 1. ∂τ l(τ ) dτ ∂ξ ∂ξ ∂τ
(32)
M(ξ, 0) = 1, T (ξ, 0) = 0, σ (ξ, 0) = 1,
(33)
∂T ∂M = 0 and = 0, ξ = 0, ∂ξ ∂ξ
(34)
∂σ ∂M ∂T −χ +λ = l(τ )Bi(T sur − 1), ξ = 1. (35) ∂ξ ∂ξ ∂ξ The further condition is needed on the moving interface to determine the position of the interface itself (see [5]). In this model it is M = 0 and −
dl M0 ∂ M = , ξ = 1, (36) dτ l(τ ) ∂ξ with l(0) = 1. Eq. (36) is known as the ‘Stefan condition’ and it expresses the heat balance on the interface. Eqs. (30)–(32) were solved together with (36), which involves a shrinkage equation obtained using the IMEX method as above. Discretizing Eq. (30)–(32) with a central difference approximation gives, ( n+1 n+1 n+1 ) n+1 n M i−1 − 2M i + M i+1 Mi − Mi 1 = 2 ∆τ l (n + 1) (∆ξ )2 ( )( n+1 n+1 ) M i+1 − M i−1 ξi l(n + 1) − l(n) n + (37) − k2 T i , l(n + 1) ∆τ 2∆ξ ( n+1 n+1 n+1 ) n T i−1 − 2T i + T i+1 − Ti Le = 2 ∆τ l (n + 1) (∆ξ )2 ( n n ) σ i−1 − 2σ in + σ i+1 Dvnon + 2 l (n + 1) (∆ξ )2 ( )( n+1 n+1 ) T i+1 − T i−1 l(n + 1) − l(n) ξi , + l(n + 1) ∆τ 2∆ξ
n+1
Ti
(38)
( )( n+1 n+1 n ) σ n − σ i−1 T − T i−1 σ in+1 − σ in ξi l(n + 1) − l(n) = µnon i+1 + i+1 ∆τ l(n + 1) ∆τ 2∆ξ 2∆ξ n+1
−
E non σin
− µnon
Ti
n
− Ti , ∆τ
(39)
respectively, for i = 1, . . . , N − 1 and n = 0, . . . , N − 1. Discretization of Eq. (36) by using the method proposed by Furzeland in [5] we have ( ) n n n 3M N − 4M N −1 + M N −2 l(n + 1) − l(n) M0 = , (40) ∆τ l(n + 1) 2∆ξ for n = 0, . . . , N − 1.
260
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265 Table 2 Input parameters used in the simulations of drying. Parameter
Value
D α Dv E µ∗ Hv Tair k1 k
8 × 10−10
m2 /s
4 × 10−9 m2 /s 8 × 10−12 (m2 ◦ C)/(Pa s) 1 × 10−9 Pa 1 × 10−4 Pa s 2345 kJ/kg 60 ◦ C 1 × 10−10 (w.b)/(s ◦ C) 0.475 W/m ◦ C
Parameter
Value
h M0 T0 σ0 L0 ∆x ∆t A k¯
30 W/m2 ◦ C 0.8 w.b 0◦ C 5 × 10−2 Pa 5 × 10−3 m 2 × 10−2 2 × 10−7 1 × 10−4 (1/◦ C) 1(1/◦ C)
Table 3 Errors and convergence orders for the temperature and moisture for different step-size ∆x. ∆x
E∆x (T )
R(T )
E∆x (M)
R(M)
0.01 0.005 0.004 0.002
3.2 × 10−3 1.6 × 10−3 1.3 × 10−3 5.1303 × 10−4
– 1 0.9305 1.3414
1.57 × 10−2 1.76 × 10−2 1.64 × 10−2 9.7 × 10−3
– −0.1648 0.3165 0.7576
Table 4 Errors and convergence orders for the temperature and moisture for different step-size ∆t. ∆t
E∆t (T )
2 × 10−3
3.68 × 10−2
2 × 10−4 2 × 10−5 2 × 10−6
8.2 × 10−3 1 × 10−3 9.9433 × 10−5
R(T )
E∆t (M)
R(M)
– 0.6520 0.9138 1.0025
2.9629 × 10−6
– 0.7874 0.9143 1.0274
4.8337 × 10−7 5.8877 × 10−8 5.5276 × 10−9
4. Numerical simulations In this section we obtain some numerical simulations for Eqs. (1)–(6) without shrinkage effect and with shrinkage effect by using Eqs. (24)–(29) and (37)–(40), respectively. For drying simulation, the input parameters are given in Table 2 for generic drying conditions for fruit [12] and parameters which are specially used for viscoelastic behavior of tomato. We start by analyzing the numerical convergence properties of the method (37)–(40). The behavior of the error E∆x and convergence order R are illustrated for different step-size ∆x when the time step-size is fixed at time t f = 1 × 10−3 . The error was computed using a reference solution obtained with ∆x = 0.001 and ∆t = 2 × 10−7 for T and M. In Table 3, E∆x is defined by ⎛ ⎞1 2 ∑ n n n n 2⎠ ⎝ E∆x (c) = ∥c∆x − c¯∆x ∥ = (c∆x ( p) − c¯∆x ( p)) , (41) p∈Ω
where
n c¯∆x
is the reference solution, and the rate R is defined by
R(c) =
E∆x1 (c) ) ∆x2 (c)
ln( E
1 ln( ∆x ) ∆x
,
(42)
2
where ∆x1 and ∆x2 represent two consecutive step sizes. We note that the numerical results presented in Table 3, suggest the convergence of the method (37)–(40) with first order convergence rate in space. Table 4 illustrates the behavior of the error E∆t and convergence order R, replacing ∆x by ∆t in Eqs. (41)– (42), for different time step-sizes when the step-size ∆x is fixed at time t f = 0.5. The error was computed using a reference solution obtained with ∆x = 0.02 and ∆t = 2 × 10−7 for T and M.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
261
Fig. 1. Profile of moisture through the sample with elapsed time τ = 0 − 1 in step ∆t = 2 × 10−7 for (a) Non-shrinkage model (b) shrinkage model (c) location of the surface.
Fig. 1(a, b) shows the moisture profile through the sample of the tomato with increasing time for the shrinkage and non-shrinkage models. Moisture decreased as time increased, but this process was a little slower for the nonshrinkage model compared to the shrinkage model. Fig. 1(c) shows the location of the free surface during shrinkage; the shrinkage of the product is shown by the change in the thickness of the tomato. The final value of shrinkage in this model can be determined readily as only solid mass is left at the end of drying, which in this case is 0.4576. Fig. 2 shows a temperature profile in the sample as a function of both position and time during drying for the non-shrinkage and the shrinkage model. The temperature at each location increases with the drying time, which is due to the higher drying temperature. Due to the difference between the air temperature and the sample temperature, the temperature profile rises rapidly in the early period of heating. As the heating period progresses, the rise in temperature attains an almost uniform profile. The overall dimension varies with time depending on the amount of moisture loss for each section, which can be seen in Fig. 2(b). The temperature distribution increased much faster for the shrinkage model compared with the non-shrinkage model. Fig. 3 represents the comparison of moisture and temperature gradients at the center of the sample, with and without the shrinkage effect. It can be seen that the temperature with shrinkage predicts a higher temperature at a specific time, compared to that with non-shrinkage. However, the predicted moisture content decreased faster when the shrinkage effect was included in the model. As a result of the thickness of the sample and decreases due to shrinkage, the moisture has less distance to cover and hence it reaches the surface faster, with the result that the time taken to dry a sample becomes shorter. Therefore, the model with shrinkage needs a shorter time for drying. Fig. 4 shows the profile of moisture in the sample at time τ = 0.25 in step of 2 × 10−7 , when different air temperatures are used during drying. An increase in drying temperature can influence the drying kinetics that can cause a decrease in the moisture of sample. We plot, in Fig. 5, the dependence of temperature at the center on the Dvnon which is proportional to the viscoelastic diffusion coefficient, Dv . In this figure we observe that the sample is acting as a barrier to the entrance
262
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Fig. 2. Profile of temperature through the sample with elapsed time τ = 0 − 1 in step 2 × 10−7 for (a) Non-shrinkage model (b) shrinkage model.
Fig. 3. Comparison of moisture and temperature ratio at the center for shrinkage and non-shrinkage model.
Fig. 4. Comparison of moisture at time τ = 0.25 in step of 2 × 10−7 with different air temperatures (for shrinkage model).
of heat. This is because of the viscoelastic properties of the tomato. According to this description, an increase in Dv the temperature at the center of the sample is decreasing and it causes the moisture not decreasing. The detection of this region is a very important aspect in terms of food safety. Fig. 6 shows that by increasing the value of Le, the temperature of the center increases. This is consistent with the result of [12], that is, for the value of Le more than 100, the interior temperature is approximately equal to the air temperature. Fig. 7 shows that by increasing the value of Bi, the temperature increases rapidly.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
263
Fig. 5. Influence of the Dv , viscoelastic diffusion coefficient, on the temperature at the center of the sample.
Fig. 6. Temperature at the center for different values of Le.
Fig. 7. Temperature at the surface for different values of Bi.
Fig. 8 shows the experiment measurement of shrinkage L/L 0 taken from experiment drying of tomato [6] at 60 ◦ C, which was compared with the simulation model of shrinkage. The thickness is assumed to decrease linearly
264
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
Fig. 8. Comparison between experiment data with simulation study of variation of shrinkage with respect to the moisture content.
with a reduction in moisture content. The trend deviates with higher moisture content but there is good agreement at lower moisture content. 5. Conclusion A model to simulate the drying of tomato has been studied. The whole process is described by a set of partial differential equations that take into account the transfer of heat from air towards the product center, diffusion of moisture outwards towards the surface and stress driven diffusion. The stability of the model subject to nonshrinkage was established under a mathematical relation between the parameters involved in the equations. To solve the system of equations an IMEX method was proposed and its convergence was numerically studied. The numerical simulations show qualitative agreement with the physical expected behavior. Moisture decreases more quickly in the shrinkage model, which concludes that the model with shrinkage effect needs a shorter time for drying. The viscoelastic properties of the tomato are shown to be an effective control mechanism to delay the heat transfer. In future work we are planning to complete the model by considering three dimensionality and studying the stability analysis in time-varying domain. From an analytical point of view the convergence analysis of the IMEX methods used in the simulations will deserve our attention in the future. References [1] C.T. Akanbi, R.S. Adeyemi, A. Ojo, Drying characteristics and sorption isotherm of tomato slices, J. Food Eng. 73 (2006) 157–163, http://dx.doi.org/10.1016/j.jfoodeng.2005.01.015. [2] D. Arslan, M.M. Ozcan, Drying of tomato slices: changes in drying kinetics, mineral contents, antioxidant activity and color parameters, CytA-J. Food 9 (2011) 229–236, http://dx.doi.org/10.1080/19476337.2010.522734. [3] E. Azhdari, J.A. Ferreira, P. de Oliveira, P.M. da Silva, Diffusion, viscoelasticity and erosion: Analytical study and medical applications, J. Comput. Appl. Math. 275 (2015) 489–501, http://dx.doi.org/10.1016/j.cam.2014.01.025. [4] B. Broyart, G. Trystram, Modelling of heat and mass transfer phenomena and quality changes during continuous biscuit baking using both deductive and inductive (neural network modelling principles), Food Bioprod. Process. 81 (2003) 316–326, http://dx.doi.org/10. 1205/096030803322756402. [5] J. Crank, Free and Moving Boundary Problem, New York Clarendon-Oxford Science Publication, 1984. [6] B. Dianda, M. Ousmane, T.K. S. Kam, D.J. Bathiébo, Experimental study of the kinetics and shrinkage of tomato slices in convective drying, Afr. J. Food Sci. 9(5) (2015) 262–271, http://dx.doi.org/10.5897/AJFS2015.1298. [7] S.K. Ghosh, G.C. Shit, J.C. Misra, Heat transfer in hydromagnetic fluid flow: Study of temperature dependence of fluid viscosity, J. Appl. Fluid Mech. 7 (2014) 633–640. [8] T. Gronwall, Note on the derivative with respect to a parameter of the solutions of a system of differential equations, Ann. of Math. 20 (1919) 292–296, http://dx.doi.org/10.2307/1967124. [9] R. Jakman, A.G. Marangoni, D.W. Stanley, The effects of turgor pressure on puncture and viscoelastic properties of tomato tissue, J. Text. Stud. 23(4): (1992) 491–505, http://dx.doi.org/10.1111/j.1745-4603.1992.tb00036.x. [10] R.G. Owens, T.N. Phillips, Computational Rheology, Imperial College Press, 2002.
E. Azhdari and A. Emami / Mathematics and Computers in Simulation 166 (2019) 253–265
265
[11] N.A. Shahari, Mathematical Modelling of Drying Food Products: Application to Tropical Fruits (Ph.D. Thesis), University of Nottingham, 2012. [12] N.A. Shahari, S. Hibberd, Analysis of single phase moisture and heat model of food drying, in: 4th International conference of Mathematical Models in Engineering and Computer science, 2013, pp. 137–142. [13] V. Thomée, B. Wahlbin, Long-time numerical solution of a parabolic equation with memory, Math. Comp. 62 (1994) 477–496. [14] I. Verlent, M. Hendickx, P. Rovere, P. Moldenaers, A.V. Loey, Rheological properties of tomato based products after thermal and high-pressure treatment, J. Food Sci. 71(3) (2006) 5243–5248, http://dx.doi.org/10.1111/j.1365-2621.2006.tb15648.x. [15] G. Xanthopoulos, S. Yanniotis, E. Talaiporou, Influence of salting on drying kinetics and water diffusivity of tomato halves, Int. J. Food Prop. 15(4) (2012) 847–863, http://dx.doi.org/10.1080/10942912.2010.506018.