transportul umiditatii

Upload
vele2ionut 
Category
Documents

view
228 
download
0
Transcript of transportul umiditatii
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 1/20
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 2/20
from air to the dough surface in the absence of baking tins, conduction in the continuousliquid/solid phase of the dough, and evaporationcondensation in the gas phase of the dough.
De Vries et al. [9] described a 4step mechanism for the heat transport inside dough: (1)water evaporates at the warmer side of a gas cell that absorbs latent heat of evaporation;(2) water vapour then migrates through the gas phase; (3) when meeting the cooler sideof the gas cell, water vapour condenses and becomes liquid; (4) finally heat and water aretransported by conduction and diffusion through the gluten gel to the warmer side of thenext cell. The water diffusion mechanism becomes important to heat transfer, because doughtends to be a poor conductor that limits the heat transfer rate via conduction.
The above described mechanism of diffusion together with evaporation and condensationin dough was subsequently adopted by Tong and Lund [22] and Thorvaldsson and Janestad [21]. With this mechanism, liquid water moved towards the loaf centre as well as thesurface by evaporation and condensation, reducing the partial water vapour pressure due to
the temperature gradient. As a result, crumb temperature change was accelerated.In the model by Zanoni et al. [23], an evaporation front inside the dough was assumed
to be always at 100oC. This evaporation front progressively advanced towards the centre asbread’s temperature increased. Crust was formed in the bread portion above the evaporationfront. With similar parameters, Zanoni et al. [24] developed a 2dimensional axisymmetricheat diffusion model. The phenomena were described separately for the upper and lowerparts (crust and crumb). The upper part (crust) temperature was determined by equationsincluding heat supply by convection, conductive heat transfer towards the inside, and convective mass transfer towards the outside. The lower part (crumb) temperature was determinedby the Fourier’s law. In addition to the Cartesian coordinate models, a 1dimensionalcylindrical coordinate model was also established by De Vries et al. [9].
Among the various models, the internal evaporationcondensation mechanism well explains the fact that heat transfer in bread during baking is much faster than that describedby the conduction alone in dough/bread. It also supports the observation that there isan increase in the liquid water content at the centre of bread during the early stage of baking (Thorvaldsson and Skjoldebrand [20]), rather than a monotonous decrease resultedfrom having liquid water diffusion and surface evaporation only. Therefore, a promisinglygood model for bread baking might be a multiphase model which consists of three partialdifferential equations for the simultaneous heat transfer, liquid water diffusion and watervapour diffusion respectively, together with two algebraic equations describing water evaporation and condensation in the gas cells. Indeed, Thorvaldsson and Janestad [21] used this
multiphase model to describe a onedimensional case where baking tin was absent.In [21], a slab of bread crumb being baked in a conventional oven was considered. The
model assumed that the vapor and liquid water diffused separately and the phase change(evaporation and condensation) occurred instantaneously, i.e., the vapor content V wasdirectly proportional to the saturated vapor content V s at the same temperature at anygiven time, provided that there was enough liquid water available. The authors showedreasonably good agreement between the numerical results predicted by the model and theexperimental measurements. However, further investigations using the same model revealed
2
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 3/20
that the numerical procedure for solving the model became a serious issue as the time stepsize was refined [18, 25].
Zhou [25] demonstrated that numerically solving the model presented a big challenge.Although various schemes of finite difference methods (FDM) and finite element methods(FEM) were applied, the corresponding solutions were all shown to be highly sensitive tothe time step and satisfactory results were only yielded over a limited range of time step.Erroneous and/or divergent results were produced when the time step size was either toosmall or too large. While it is reasonable to expect that large time step may lead to instability,failure of the numerical methods due to small time step sizes is nonintuitive.
The main objective of this paper is to provide an explanation as well as a remedy tothis highly nonintuitive outcome of the instantaneous phase change model proposed in [21].In order to identify the source of instability in the instantaneous phase change model, weconstruct a reactiondiffusion model by assuming that the phase change follows a modified
HertzKnudsen equation, which can be derived using statistical mechanics principles [13].We show that the two models are related from a numerical point of view, when a specialtime stepping scheme is used for the reactiondiffusion model. The instantaneous phasechange model is equivalent to the reactiondiffusion model with a variable evaporation ratereciprocal to the time step size. Therefore, reducing the time step size in the instantaneousphase change model has the additional effect of refining the time step as well as increasingthe evaporation rate in the reactiondiffusion model.
The reactiondiffusion model allows us to study the effects of the evaporation rate in moredetails by separating two types of instabilities, i.e., instability associated with the particularnumerical scheme and an intrinsic diffusive instability associated with the model. Combining numerical tests and linear stability analysis, we show that the diffusive instability is anintrinsic feature of the model, which does not disappear as the time step is reduced. As aresult, the twophase region where vapor and liquid water coexists, may become unstable.For the parameter values used in [21, 18, 25], the constant state solution in the twophaseregion is linearly unstable even though the solution is stable when the phase change anddiffusion are considered separately. Furthermore, the rate of growth depends on the evaporation rate. A larger evaporation rate leads to faster growth of the disturbances. Oscillationin solutions can occur in the twophase region before the dry region completely takes over.
We show that the reactiondiffusion model does not lead to numerical instability whensufficiently small time step size is used. In the instantaneous phase change model, reducingthe time step size is equivalent to increasing the evaporation rate while keeping the product
of the time step size and the evaporation rate fixed. Therefore, systematic refinement of thetime step size induces faster growth of error, which eventually causes numerical instabilityobserved in the computations in [18, 25].
In order to put the problem we study in this paper in a broader context, we shall givea brief overview of some relevant applications whose mathematical description has certainsimilarity to the bread baking models discussed here. A general feature of these problems isthe coupling of thermal diffusion and phase change. We refer to [3, 5, 6, 7, 10] for treatmentsof phase change and heat transfer phenomena in general. A closely related application is
3
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 4/20
given in [16, 15], where a moisture transport model for the wetting and cooking of a cerealgrain is considered. The temperature is decoupled from the moisture model and it is used as aparameter in their model. In [12, 11], the condensed phase combustion or gasless combustionis considered. The model has applications in synthesizing certain ceramics and metallic alloysand involves a reactiondiffusion system for the temperature and the concentration of thefuel, where the diffusion coefficient of the fuel concentration is set to be zero. In [2], amodel for the aggregate alkali reaction in fluid leaching processes, which is similar to the thegasless combustion, is studied numerically. All these models are similar to the one studiedin [21, 18, 25], with some differences in the way the phase change is handled and in theparameters and coefficients. Numerical techniques used in those studies are also similar,where either finite difference and finite element methods or pseudospectral method are usedfor spatial discretization.
The rest of the paper is organized as follows. In section 2, we propose a reactiondiffusion
model based on the HertzKnudsen equation. The relationship between our model and theinstantaneous phase change model is explored in section 2.2 when the numerical method isdiscussed. Numerical results for the reactiondiffusion model are given in section 2.3. Insection 3 we carry out the linear stability analysis of our new model. We finish the paper bya conclusion and a short discussion on future directions in section 4.
2 A ReactionDiffusion Model
Following [21, 18, 25] we assume that the bread slab can be treated as a one dimensionalhomogenous porous medium (0 < x < L) with density ρ, specific heat capacity c p, thermal
conductivity k. While in reality both c p and k depend on water content, they are assumedto be constants in [21], where density ρ is given as a linear function of the liquid watercontent. The main variables are temperature T , liquid water content W and water vaporcontent V (defined as the percentage of liquid water and vapor mass, compared to that of the mixture including the bread) inside the bread with respect to the total weight. Thus thetotal vapor and liquid water masses per unit of volume can be computed as ρV and ρW ,respectively. The vapor and liquid water can be generated via evaporation and condensation,respectively, determined by the saturated vapor concentration V s (or saturation pressure P s),which is temperature dependent. The vapor and liquid water can also be transported viadiffusion, with coefficients Dv and Dw, respectively.
The governing equations for T , V and W are
ρc p∂T
∂t=
∂
∂x
k
∂T
∂x
+ λΓ, (1)
∂V
∂t=
∂
∂x
Dv
∂V
∂x
− Γ
ρ, (2)
∂W
∂t=
∂
∂x
Dw
∂W
∂x
+
Γ
ρ. (3)
4
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 5/20
Here Γ is the rate of the phase change (mass per unit volume per unit time), given bythe modified HertzKnudsen equation (See [13] where this equation is used for studyingevaporation phenomenon)
Γ = E (1− φ)
M
2πR
(P v − cP s)√T
(4)
where E is the condensation/evaporation rate, φ is the porosity of the bread slab, M isthe molecular weight of water and R is the universal gas constant, P v and P s are the vaporpressure and saturation pressure, and c is another phase change constant same as the modelin the previous section. For simplicity we have included the pore surface area in E , which isinversely proportional to the pore size of the porous sample. Even for a unit pore surface,the value of evaporation rate E is subject to debate [13]. In this study, we will estimate the
value of E by comparing the current model with the instantaneous phase change one in thenext section.Assuming the ideal gas law,
P v =RρV T
φM , P s =
RρV sT
φM (5)
we can rewrite (4) as
Γ =E (1− φ)ρ
φ
RT
2πM (V − cV s) . (6)
It is worth noting that evaporation can only occur when there exists sufficient amount of
liquid water. Thus, if V < cV s and W = 0, then no evaporation occurs and Γ = 0.As in [21, 18, 25], we assum that the bread has certain moisture content W 0 (liquid water)
at the room temperature T 0 initially. At time zero, the bread is placed in an oven preheatedat temperature T a by a radiator at temperature T r. The boundary conditions are mixedcondition at x = 0:
−k∂T
∂x= hr(T r − T ) + hc(T a − T ), (7a)
−∂V
∂x= hv(V a − V ), (7b)
−∂W
∂x
= hw(W a−
W ) (7c)
and symmetric condition at x = L:
∂T
∂x=
∂V
∂x=
∂W
∂x= 0. (7d)
Here hr, hc, hv and hw are the radiative, convective heat transfer coefficients, vapor andliquid water mass transfer coefficients. V a and W a are the vapor and liquid water content inthe oven air.
5
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 6/20
2.1 Nondimensionalization
We now proceed to nondimensionalize the equations by choosing the following scaling
θ =T − T 0
T , τ =
t
t, ξ =
x
L, ρ′ =
ρ
ρ.
Here T 0 = 298.15 K is the initial temperature and T = T r − T 0 is the difference between theradiator and initial temperature. ρ = 284 kg/m3 is the density of the flour. We choose thediffusive time scale for the temperature equation
t =ρc pL2
k.
Drop the prime for simplicity and we have the following equations
θτ =1
ρθξξ + S, (8)
V τ = D1[(θ + θ0)2V ξ]ξ − αS, (9)
W τ = D2W ξξ + αS (10)
where the nondimensional density is given as ρ = W + ρ0. The other dimensionless parameters are defined as
D1 =DV 0ρT 2c p
k, D2 =
DW ρc pk
, α =T c p
λ.
When there is sufficient amount of liquid water, the source term for the phase change isgiven byS = β
θ + θ0 (V − cV s) (11)
where
β = λE (1− φ)ρL2
φk
R
2πM T ,
and the saturated vapor content is
V s =V s,0eγθ − V s,1
(ρ0 + W )(θ + θ0),
with parameters
V s,0 =P s,0M
ρT R, V s,1 =
P s,1M
ρT R, γ = κT .
Note that we have fitted the saturated vapor pressure data [19] by an exponential function
P s = P s,0 exp[κ(T − T 0)]− P s,1. (12)
6
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 7/20
The boundary conditions at x = 0 (ξ = 0) and x = L (ξ = 1) are nondimensionalizedsimilarly. At ξ = 0, we have
θξ = (h3 + h4)(θ − 1), (13a)V ξ = h1(V − V a), (13b)
W ξ = h2(W −W a) (13c)
where
h1 = hV L,
h2 = hW L,
h3 =σT 3L
k
[(1 + θ0)2 + (θ + θ0)2](θ + 1 + 2θ0)
ǫ−1 p + ǫ−1
r − 2 + F −1sp
,
h4 = hcLk
.
Here we have assumed that the air temperature T a equals to the radiator temperature T r.At ξ = 1, we have
θξ = V ξ = W ξ = 0. (13d)
2.2 Numerical Method
We now turn our attention to the numerical procedure by describing a splitting scheme forthe reactiondiffusion model (cf. [4]). We will revisit the instantaneous phase change model
and establish a connection between our model and the instantaneous phase change model.Numerical solutions of our new model will be presented at the end of the section.
2.2.1 Reactiondiffusion model
Splitting scheme. We adopt a splitting method for the reactiondiffusion model by separating the diffusion from the phase change (reaction) term and solving the equations in twosteps 1. We first solve the vapor and liquid water equations by
V τ = −αS, (14)
W τ = αS. (15)
The second equation for liquid water can be replaced by an algebraic constraint
V + W = V 0 + W 0
1From the numerical method point of view, there is also a benefit to use the splitting. Since we canobtain an explicit relation between W and V in the phase change stage, W can be computed using V
without treating the source term S in the W equation (or with that relation the source term in the W
equation can be written in terms of W and is dissipative). As a result, the algorithm is more stable. Wenote that time splitting schemes have been used in the literature for reactiondiffusion equations and we referinterested readers to [1, 8] and references therein.
7
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 8/20
where V 0 and W 0 are the values at the beginning of the splitting step. Diffusion of vaporand liquid water as well as the temperature equation are solved in the second step by
θτ = 1ρ0 + W
θξξ + S, (16)
V τ = D1[(θ + θ0)2V ξ]ξ, (17)
W τ = D2W ξξ . (18)
Time stepping scheme. The numerical procedure in semidiscrete form can be describedas follows.
1. We solve for vapor content first using
V c
−V (n)
∆τ =−
αS 1
(19a)
whereS 1 = β
θ(n) + θ0(V − cV (n)s ). (19b)
Here superscript c denotes the solution updated due to phase change alone and V isthe arithmetic average, i.e.,
V =V c + V (n)
2. (20)
2. Evaporation can only occur if there exists sufficient amount of the liquid water. In otherwords, liquid water content must remain nonnegative. If the available water W (n) is
less than the amount V c
−V (n)
computed by (19), then the liquid water becomes zero.Otherwise, the amount of the liquid water is given by W (n) − V c + V (n). Therefore,
W c = max{W (n) + V (n) − V c, 0}. (21)
3. As a consequence, vapor content needs to be corrected using the constraint
V c = W (n) + V (n) −W c. (22)
4. To account for the possibility of running out of liquid water, phase change needs to becorrected using
S 2 = W c
−W (n)
α∆τ . (23)
5. V (n+1) and W (n+1) are updated using the diffusion equations
V (n+1) − V c
∆τ = (DvV
(n+1)ξ )ξ, (24)
W (n+1) −W c
∆τ = (DwW
(n+1)ξ )ξ. (25)
8
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 9/20
6. Update the temperature using the rate of phase change S 2 given by (23)
θ(n+1)
−θ(n)
∆τ =
1
ρ0 + W (n+1) θ
(n+1)
ξξ + S 2. (26)
The time discretizations described above are standard and it is straight forward to verifythat the scheme is consistent and formally the order of accuracy is ∆ τ for reasonably chosenV , such as the arithmetic average used in this paper.
2.2.2 Instantaneous phase change model in [21]
The instantaneous phase change model used and studied in [21, 18, 25] is essentially a discretetime model which can be described as follows (in the nondimensional form).
1. We update vapor and water contents using
V ∗ = min{cV (n)s , V (n) + W (n)}, W ∗ = V (n) + W (n) − V ∗ (27a)
orW ∗ = max{W (n) + V (n) − cV (n)s , 0}, V ∗ = V (n) + W (n) −W ∗. (27b)
2. Rate of phase change is computed using
S 1 =W ∗ −W (n)
α∆τ . (28)
3. Update vapor due to diffusion
V ∗∗
−V ∗
∆τ = (DvV ∗∗
ξ )ξ. (29)
4. We update vapor and water contents again using
V (n+1) = min{cV (n)s , V ∗∗ + W ∗}, W ∗∗ = V ∗∗ + W ∗ − V (n+1) (30a)
orW ∗∗ = max{W ∗ + V ∗∗ − cV (n)s , 0}, V (n+1) = V ∗∗ + W ∗ −W ∗∗. (30b)
5. Rate of additional phase change is computed using
S 2 =W ∗∗ −W ∗
α∆τ . (31)
6. W (n+1) is updated using the diffusion equations
W (n+1) −W ∗∗
∆τ = (DwW
(n+1)ξ )ξ. (32)
7. Update the temperature
θ(n+1) − θ(n)
∆τ =
1
ρ0 + W (n+1)θ(n+1)ξξ + S 1 + S 2. (33)
9
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 10/20
2.2.3 Simplified instantaneous phase change model
The rate of phase change given by (31) is typically small compared to that from (28).
Therefore, we can simplify the model by eliminating steps (30) and (31) and using thefollowing procedure.
1. We update vapor and water contents using
W ∗ = max{W (n) + V (n) − cV (n)s , 0}, V ∗ = V (n) + W (n) −W ∗. (34)
2. Rate of phase change is computed using
S =W ∗ −W (n)
α∆τ . (35)
3. Update vapor due to diffusion
V (n+1) − V ∗
∆τ = (DvV
(n+1)ξ )ξ. (36)
4. W (n+1) is updated using the diffusion equations
W (n+1) −W ∗
∆τ = (DwW
(n+1)ξ )ξ. (37)
5. Update the temperature
θ(n+1) − θ(n)∆τ
= 1ρ0 + W (n+1)
θ(n+1)ξξ + S. (38)
2.2.4 Model comparison
Note that the diffusion part of the reactiondiffusion model (24)(26) are the same as thatfor the simplified instantaneous model (36)(38). The only difference between the reactiondiffusion model and the instantaneous phase change model lies in the way the phase changeis computed, or more precisely the way vapor content is computed. We now show that thetwo models are related to each other in the following sense.
In the instantaneous phase change model, assuming there is sufficient amount of water 2,
we setV ∗ = cV s. (39)
Again assuming there is sufficient amount of water, using the reactiondiffusion model andV defined in (20), we obtain
V c − V (n) = −ν (V c + V (n) − 2cV (n)s )
2When there is not sufficient amount of water, it is straightforward to show that the two models areequivalent since V
c = V ∗ = V (n) + W
(n) and W c = W ∗ = 0.
10
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 11/20
where
ν =∆τ αβ
√θ + θ0
2.
When ν = 1, or
β =2
∆τ α√
θ + θ0, (40)
which yieldsV c = cV (n)s .
Therefore, we show that the instantaneous phase change model is numerically equivalentto the reactiondiffusion model with a variable β inversely proportional to the time stepsize ∆τ . From numerical analysis point of view, this suggests that the instantaneous phasechange model can be viewed as an inconsistent discretization of the reactiondiffusion model,
which provides a possible explanation for the observed divergence in the numerical solutionas ∆t → 0.
2.3 Numerical Results
The numerical solutions are obtained using parameter values given in [21, 25]. To simulatethe moisture transport during bread baking we need to determine the evaporation/condensationrate E or its nondimensionalized quantity β . This is done by choosing its value so that thenumerical solution matches the experimental results in [21]. The estimated value is β = 781and the corresponding evaporation rate is E = 4.5 × 10−3. In Figure 1, the numerical results based on this evaporation rate is given for two time step sizes ∆t = 10 and 0.2 second
using a coarse grid (16 grid points in x) and a fine grid (128 grid points in x), respectively.The numerical procedure remains stable as the time step is refined for a fixed grid size inx, contrary to the instantaneous phase change model. We note that the nonsmoothnessof the solution in Figure 1a is due to the coarse spatial and temporary grids used in thecomputation. When we refined the spatial grid, the solution becomes smooth, as shown inFigure 1b.
Recall from (40) that reducing the time step size in the instantaneous phase change modelis equivalent to increase the value of E in the reactiondiffusion model. Therefore, it willhelps us to understand the mechanism of the instability associated with the instantaneousphase change model by carrying out the computations using a larger value of E . In Figure 2,computational results using the reactiondiffusion model with E = 0.0045 and E = 0.045are presented. Here we plot the snap shots of the water content at various times in orderto show what can happen with the model if we increase E . The results are obtained using128 grid points in the x variable and a time step size of 0.2 second. We have also donefurther refinement tests but the results remain virtually the same, indicating convergence of the numerical solutions. Since we have decoupled the evaporation rate from the time stepsize, we can use a much smaller time step size while keeping the evaporation rate unchanged,which is not possible for the instantaneous phase change model. From Figure 2, it can be
11
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 12/20
0 20 40 60 80 1000
0.05
0.1
V
0 20 40 60 80 1000
0.5
1
W
0 20 40 60 80 1000
100
200
300
t (min)
T
(a)
0 20 40 60 80 1000
0.05
0.1
V
0 20 40 60 80 1000
0.5
1
W
0 20 40 60 80 1000
100
200
300
t (min)
T
(b)
Figure 1: Computational results of V , W , and T obtained by using the reactiondiffusion model: (a). Time and spatial step sizes are ∆t = 10 second and ∆x = 2−4 × 10−2 cm; (b).Time and spatial step sizes ∆t = 0.1 second and ∆x = 2−8× 10−2 cm. Here the dotted linesrepresent the quantities at the surface of the slab ( x = 0), the dashed lines are for thoselocated in the middle of the domain ( x = L/2) and the solid lines represent the values at thecenter of the slab ( x = L). The time is measured in minutes in order to compare our resultswith those from previous studies in [21, 18, 25].
12
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 13/20
0
0.002
0.004
0.006
0.008
0.01
0
10
20
30
40
50
60
70
80
0
1
2
t ( mi n )
x (m)
W
0
0.002
0.004
0.006
0.008
0.01
0
5
10
15
20
25
30
35
40
0
1
2
t ( mi n )
x (m)
W
(a) (b)
Figure 2: Snap shots of the water content computed using the reactiondiffusion model for two values of the evaporation constant: (a). E =0.0045; (b). E =0.045. For the smaller value of E , the water content remains monotonic in x and the interface between the dry and wet regions is clearly defined. For the larger value of E , disturbance from the interface growsand eventually leads to multiple drywet regions.
seen that the solution for the larger evaporation rate starts to oscillate. However, this is notdue to numerical instability.
We have also experimented with the instantaneous phase change model and the resultsconfirmed the observations in [18, 25]. Due to space limitation, we will not present thoseresults here and instead refer interested readers to [18, 25] for more details. In the next sectionwe will provide an explanation for the observed oscillation in the solution of the reactiondiffusion model associated with relatively large value of E (or small ∆t in the instantaneousphase change model when the numerical procedure is still stable).
3 Linear Stability Analysis
We now turn our attention to the stability of the solution of the nondimensional model (8)
(10) near a steady state in an infinite domain. It is easy to see that the constant state V 0,W 0 and θ = 0 satisfies the equations as long as
V 0 =c (V s,0 − V s,1)
θ0(ρ0 + W 0).
To examine the stability of this constant state solution, we carry out linear stabilityanalysis by assuming that
θ = θ exp(sτ + imξ), V = V 0 + v exp(sτ + imξ), W = W 0 + w exp(sτ + imξ).
13
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 14/20
The equations for θ, v and w are
sθ = −m2
ρ0 + W 0 θ + β −V 1θ + v +
V 0ρ0 + W 0 w
, (41)
sv = −θ20D1m2v − αβ
−V 1θ + v +
V 0ρ0 + W 0
w
, (42)
sw = −D2m2w + αβ
−V 1θ + v +
V 0ρ0 + W 0
w
, (43)
where β =√
θ0β and
V 1 =c (V s,0 − V s,1)
θ0(ρ0 + W 0)− cV s,0γ
θ0(ρ0 + W 0).
Rearranging (41)(43) in the matrix form, we obtain
M y = sy
where
y =
θ
vw
, M =
− m2
ρ0+W 0− βV 1 β β V 0
ρ0+W 0
αβV 1 −θ20D1m2 − αβ −αβ V 0ρ0+W 0
−αβV 1 αβ −D2m2 + αβ V 0ρ0+W 0
. (44)
3.1 A Special Case
We want to find out whether instability would occur in this system. This requires us to findeigenvalues of (44). To avoid complicated calculations we consider a special case where thelinearized diffusion coefficients for temperature and liquid water are same, i.e.
1
ρ0 + W 0= θ20D1. (45)
For our choice of parameters in this breadbaking problem these two coefficients are of similarsizes. Therefore we expect that this assumption would not change the nature of stability of the system. Under this assumption the eigenvalues satisfy the following equation
(θ20D1m2 + s)
s2 +
D2m2 + θ20D1m2 + βV 1 + αβ − αβV 0ρ0 + W 0
s
+D2m2(θ20D1m2 + βV 1 + αβ )− αβV 0ρ0 + W 0
θ20D1m2
= 0. (46)
14
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 15/20
3.1.1 Reactiononly
When there is no diffusion, we have m = 0, and the eigenvalues are s1 = s2 = 0 and
s3 = −αβ −
V 1V 0− α
ρ0 + W 0
βV 0.
For our problem, α = 0.2864, γ = 9.662, ρ0 = 0.5986, θ0 = 1.6116 and W 0 = 0.4. Thus,s3 < 0 and the constant solution is stable. Note that in this case the assumption on thediffusion coefficients for temperature and liquid water is satisfied.
3.1.2 Reactiondiffusion
We can easily see the eigenvalue associated with the first factor of equation (46) is negative.
Note that the coefficient of s in the second factor
C 1 = D2m2 + θ20D1m2 + βV 1 + αβ − αβV 0ρ0 + W 0
> 0.
The signs of two eigenvalues associated with the second factor is determined by the sign of the coefficient of s0:
C 0 = D2m2(θ20D1m2 + βV 1 + αβ )− αβV 0ρ0 + W 0
θ20D1m2
since the eigenvalues from the second factor is given by
−C 1
± C 21
−4C 0. So if C 0 < 0 and
C 21 > 4C 0, we will have a positive eigenvalue which indicates instability of the system. If D2
is comparable to D1 then C 0 will be positive and both eigenvalues will be negative and thesystem is stable. Since in our case D2 is much smaller than D1, C 0 will be negative if m is notsufficiently large, in which case we will have one positive eigenvalue which leads to instability.
Remark. It is interesting to note that in this case diffusion is actually destablizing. Theinstability discussed here is conceptually related to Turing instability or diffusive instabilityobserved in pattern formation, crystal and tumor growth; cf. [17, 14]. However, the mathematical model and the nonlinearity related to bread baking are quite different from thoseapplications.
3.2 The General Case
In the case when diffusion coefficients for the temperature and vapor equations are not thesame, analytical expression can still be obtained. However, the formula is too complicatedto provide useful insights. We can nevertheless find the eigenvalues numerically by using aMatlab routine eig.m. Again, for the relevant parameter values two of the three eigenvalueshave nonpositive real parts. The real part of the third eigenvalue is plotted in Figure 3 whereboth the wave number and the eigenvalue are normalized as s = s/β and m = m2/β . It can
15
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 16/20
0 0.2 0.4 0.6 0.8 1−12
−10
−8
−6
−4
−2
0
2
4
6x 10
−4
Scaled Wave Number
S c a l e d G r o w t h R a t e
Figure 3: Plot of the normalized eigenvalue s v.s. the normalized wave number m. The circleindicates the most unstable wave number mmax =0.09. The square indicates the upper limit of the unstable modes, mb =0.4469.
be seen that there exists a range of wave numbers between zero and a finite value, indicatedby the square in the figure, within which the perturbation will grow. Larger wavenumberdisturbances beyond the critical wave number will not grow. Furthermore, there exists awavenumber which grows the fastest, indicated by the circle in the figure. Finally, the rangeof unstable frequencies increase linearly with
β and the growth rate increases linearly with
β . Therefore, a larger evaporation rate E (implies a larger β ) leads to a wider range of unstable modes and faster growth rates for all the unstable modes.
We now verify the linear stability analysis by solving the reactiondiffusion model withthe Dirichlet conditions which permit the constant solution θ = 0, V = V 0 and W = W 0.In Figure 4a, the computed liquid water content at the end of 100 minutes is shown. Thesolutions are obtained by superimpose a small disturbance in the form of ǫ cos[m(ξ − 0.5)]to the constant liquid water content W 0 with m = 2πω, ω =0.5, 1.5, and 4 and ǫ = 0.1. The
evaporation rate is E = 0.0045 and the range of unstable frequencies is 0 < ω < 2.9738.The most unstable mode is given by ω ≈ 1.3345 and the rate of growth is approximatelyexp (3.5056× 10−4t) measured in dimensionless time. The results show that even thoughthe disturbances with wave numbers below ω ≈ 3 grow with time, it takes a relatively longperiod of time (on the scale of 80 minutes for the fastest growing mode) to develop. Onthe other hand, the evaporation of all the liquid water takes about 60 minutes (Figures 1and 2a). Therefore, all the liquid water in the bread would have already evaporated beforethe instability develops and takes effect.
16
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 17/20
0 0.002 0.004 0.006 0.008 0.010
0.5
1
W
x (cm)
0 0.002 0.004 0.006 0.008 0.010
0.5
1
W
0 0.002 0.004 0.006 0.008 0.010
0.5
1
W
0 0.002 0.004 0.006 0.008 0.010
0.5
1
W
0 0.002 0.004 0.006 0.008 0.010
0.5
1
W
0 0.002 0.004 0.006 0.008 0.010
0.5
1
W
x (cm)
(a) (b)
Figure 4: Computational results of liquid water content subjected to disturbances. (a)E =0.0045 with ω =0.5, 1.5 and 4. (b) E =0.045 with ω =1.5, 4 and 10. The dashed lines indicate the initial states and the solid lines are the final solutions.
In Figure 4b, the computed liquid water content at the end of 10 minutes is shown fora larger evaporation rate E = 0.045, subjected to the same small disturbances. In thiscase, the range of unstable frequencies is 0 < ω < 9.4. The most unstable mode is givenby ω ≈ 4.22 and the rate of growth is approximately exp (3.5056× 10−3t) in dimensionless
time. Compared to the previous case, this leads to a more rapid growth for the most unstablemode, on the order of 10 minutes in the dimensional time, while it takes about 30 minutesfor all the liquid water to evaporate, cf. Figure 2b, where the effect of the instability isclearly visible even though the sizes of the disturbances are unknown.
The difference in the growth rates for the same disturbance under the two evaporationrates can be seen in Figure 5 where the liquid water content in the middle of the domain isplotted. The prediction using the linear stability analysis is also plotted. Note that the linearstability analysis is only valid near the constant state. Nevertheless, the two sets of data arenot too far off. Thus, subject to random perturbation, the initially constant solution in theuniform twophase region will become unstable. The fastest growing mode eventually causesa multiple region of dry and twophase region, as shown in Figure 2.
Remark. From the analysis we see that the diffusive instability occurs for any positive β (or E ). However, for relatively small β computational results in the previous section aresatisfactory. One reason is that the positive eigenvalue is small as β is relatively small andthe disturbances grow slowly. In addition, the oscillation (diffusive instability) is only visiblewhen the wave length of the unstable modes are shorter than the domain size in x. Thisexplains why no oscillation is not observed in the case of the smaller β (or E ).
17
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 18/20
0 20 40 60 80 1000
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
t (min)
W m i d ω=0.5
ω=1.5
ω
=4
0 2 4 6 8 100
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
t (min)
W m i d
ω=1.5ω=4
ω=10
(a) (b)
Figure 5: Damping rates for (a) E =0.0045 and (b) E =0.045. The liquid water content in the middle of the domain where the growth is the fastest is plotted for wave numbers ω =0.5,1.5 and 4 in (a) and ω =1.5, 4 and 10 in (b). The circles are predictions by the linear stability analysis.
4 Conclusion
The work in this paper is motivated by the simulation results in [18, 25], based on a multi
phase model for bread baking proposed in [21]. This model allows simultaneous heat, vaporand liquid water transfer by assuming that the phase change is instantaneous. However, previous study [18, 25] showed that this model only produces reasonable solutions for specificchoices of the grid and time step sizes. Reducing the grid and time step sizes usually leadto numerical instability and cause the solution to blow up.
By constructing a reactiondiffusion model and establishing a link between our modeland the model in [21], we have identified the source of the instability associated to the modelin [21]. We have shown that the instability observed in [18, 25] is a combination of two factors:the numerical instability as well as a diffusive instability. Using our reactiondiffusion model,we showed that the numerical instability can be eliminated by using a sufficiently small timestep. This is due to the fact that the reactiondiffusion model separates the numerical
instability from the diffusive instability. The diffusive instability, on the other hand, is anintrinsic feature of the model, as demonstrated by the linear stability analysis and numericaltests. For relatively large evaporation rate, the diffusive instability leads to an oscillatorysolution with multiple regions of dry and twophase zones.
Our analysis of the reactiondiffusion model also reveals that the diffusive instability isrelated to the value of the evaporation rate, which is affected by the properties of the porousmedium, such as the surface area of the pore space. This suggests that the phenomenonrelated to diffusive instability may be realized in the physical process of bread baking. Further
18
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 19/20
experimental investigation is necessary to validate or invalidate our model applied to breadbaking. On the other hand, the model itself is quite general and may be applicable to similarproblems with simultaneous heat and mass transfer processes.
Finally, we wish to point out that a consistent instantaneous phase change model maybe derived based on our reactiondiffusion model, using asymptotic analysis and β −1 as asmall parameter and this will be pursued in a future study.
Acknowledgement. We would like to thank the referees for many constructive suggestionswhich have helped to improve the quality of our manuscript. We also wish to acknowledgebeneficial discussions with Dr. Jonathan Wylie during the preparation of the manuscript.This work is partially supported by research grants from NSERC and MITACS of Canada(HH) and Academic Research Grants R146000053112 and R146000181112 from theNational University of Singapore (PL and WZ).
References
[1] D. Alemani, B. Chopard, J. Galceran, and J. Buffle, LBGK method coupled to timesplitting technique for solving reactiondiffusion processes in complex systems, Phys.Chem. Chem Phys. 7 (2005), no. 18, 33313341.
[2] G. Carey, N. Fowkes, A. Staelens and A. Pardhanani, A class of coupled nonlinearreaction diffusion models exhibiting fingering, J. Comp. Appl. Math. 166 (2004), 8799.
[3] H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford,1959.
[4] P.G. Ciarlet and J.L. Lions, Handbook of Numerical Analysis, Vol. 1, NorthHolland,Amsterdam, 1990.
[5] J. Crank, The Mathematics of Diffusion , Clarendon Press, Oxford, 1975.
[6] J. Crank, Free and Moving Boundary Problems, Oxford University Press, Oxford, 1987.
[7] H.S. Davis, Theory of Solidifications, Cambridge University Press, Cambridge, 2001.
[8] S. Descombes, Convergence of a splitting method of high order for reactiondiffusion
systems, Math. Comp. 70 (2000), no. 236, 14811501.
[9] U. De Vries, P. Sluimer and A.H. Bloksma, A quantitative model for heat transport indough and crumb during baking. In ”Cereal Science and Technology in Sweden,” ed. byN.G. Asp. STU Lund University, Lund, 1989, p.174188.
[10] A.C. Fowler, Mathematical Models in the Applied Sciences, Cambridge University Press,1997.
19
8/3/2019 transportul umiditatii
http://slidepdf.com/reader/full/transportulumiditatii 20/20
[11] M.L. Frankel, G. Kovacic, V. Roytburd and I. timofeyev, Finitedimensional dynamicalsystem modeling thermal instabilities, Physica D 137 (2000), 295315.
[12] M.L. Frankel, V. Roytburd and G. Sivashinsky, A sequence of period doublings andchaotic pulsations in a free boundary problem modeling thermal instabilities, SIAM J.Appl. Math. 54 (1994), No.4, 11011112.
[13] F. Jones, Evaporation of Water , CRC Press, 1991.
[14] D.A. Kessler, J. Koplik and H. Levine, Pattern selection in fingered growth phenomena,Adv. Phys. 37 (1988), 255339.
[15] K.A. Landman, C.P. Please, Modelling moisture uptake in a cereal grain, IMA J. Math.Buss. Ind. 10 (2000), 265287.
[16] M.J. McGuinness, C.P. Please, N. Fowkes, P. McGowan, L. Ryder and D. Forte, Modelling the wetting and cooking of a single cereal grain, IMA J. Math. Buss. Ind. 11(2000), no.1, 4970.
[17] J. D. Murray, Mathematical Biology , SpringerVerlag, Berlin, 1989.
[18] G. G. Powathil, A heat and mass transfer model for bread baking: an investigation usingnumerical schemes, M.Sc. Thesis, Department of Mathematics, National University of Singapore, 2004.
[19] R.P. Singh and D.R. Heldman, An Introduction to Food Engineering , 3rd Edition, Aca
demic Press, London 2001.
[20] K. Thorvaldsson and C. Skjoldebrand, Water diffusion in bread during baking.Lebensm.Wiss. u.Technol., 31 (1998), 658663.
[21] K. Thorvaldsson and H. Janestad, A model for simultaneous heat, water and vapordiffusion. J. Food Eng. 40 (1999) 167172.
[22] Tong, C.H. and Lund, D.B., Microwave heating of baked dough products with simultaneous heat and moisture transfer, Journal of Food Engineering , 19 (1993), 319339.
[23] B. Zanoni, C. Peri and S. Pierucci, A study of the breadbaking process I: a phenomeno
logical model, Journal of Food Engineering , 19 (1993), No. 4, 389398.
[24] B. Zanoni, C. Peri and S. Pierucci, Study of bread baking processII: Mathematicalmodelling, Journal of Food Engineering , 23 (1994), No. 3, 321336.
[25] W Zhou, Application of FDM and FEM to solving the simultaneous heat and moisture transfer inside bread during baking, International Journal of Computational Fluid Dynamics, 19 (2005), No. 1, 7377.
20