transportul umiditatii

of 20 /20
Moisture Transport and Diusive Instability during Bread Baking H. HUANG , P. LIN and W. ZHOU March 28, 2007 Abstract In this paper we study multiphase models for simultaneous heat and mass transfer process durin g bread baking. Our main objectiv e is to provid e an explanatio n and a remedy to the observed erroneous and/or divergent results associated with the instanta- neous phase change model used in the literature. We propose a reaction-diusion model based on the Hertz-Knudsen equation, where the phase change is not instantaneous but determine d by an ev aporatio n/cond ensation rate. A splitting scheme is desig ned for the reaction-diusion model so that a link between these two models can be estab- lished and the non-intuitive numerical instability associated to the instantaneous phase change model can be identied and eliminated. The evaporation/condensation rate is estimated by comparing the results of the reaction-diusion model with the experimen- tal observation s in the literature . F or the evapora tion/c ondensation rate bey ond the estimated value, oscillatory solution with multiple regions of dry and two-phase zones is observed. We show that the se are cau sed by an instabilit y intrinsic to the model (which we call diusive instability) using linear stability analysis and numerical tests. Keywords: Diusive Instability, Finite Dierence Method, Heat and Mass Transfer, Linear Stability Analysis, Multi-phase Modeling, Phase Change, Reaction-Diusion Equation. 1 Introducti on Bread baking process is dicult to model, partly due to the fact that simultaneous heat and mass transfer are involved during the process. During baking, heat transfer in dough is a combination of conduction/radiation from band or tins to the dough surface, convection Depar tment of Mathematic s and Stati stics, Y ork Univ ersity, T oronto, Ont ario, Canada M3J 1P3. [email protected] Department of Mathematics, National University of Singapore, Singapore 117543. [email protected] Food Science and Technology Programme, Department of Chemistry, National University of Singapore, Singapore 117543. [email protected] 1

Transcript of transportul umiditatii

Page 1: transportul umiditatii

8/3/2019 transportul umiditatii 1/20

Page 2: transportul umiditatii

8/3/2019 transportul umiditatii 2/20

from air to the dough surface in the absence of baking tins, conduction in the continuousliquid/solid phase of the dough, and evaporation-condensation in the gas phase of the dough.

De Vries et al. [9] described a 4-step 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 Janes-tad [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 2-dimensional axi-symmetricheat 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 convec-tive mass transfer towards the outside. The lower part (crumb) temperature was determinedby the Fourier’s law. In addition to the Cartesian co-ordinate models, a 1-dimensionalcylindrical co-ordinate model was also established by De Vries et al. [9].

Among the various models, the internal evaporation-condensation mechanism well ex-plains 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 evapo-ration and condensation in the gas cells. Indeed, Thorvaldsson and Janestad [21] used this

multiphase model to describe a one-dimensional 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


Page 3: transportul umiditatii

8/3/2019 transportul umiditatii 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 non-intuitive.

The main objective of this paper is to provide an explanation as well as a remedy tothis highly non-intuitive 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 reaction-diffusion model by assuming that the phase change follows a modified

Hertz-Knudsen 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 reaction-diffusion model. The instantaneous phasechange model is equivalent to the reaction-diffusion 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 reaction-diffusion model.

The reaction-diffusion 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. Combin-ing 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 two-phase region where vapor and liquid water co-exists, may become unstable.For the parameter values used in [21, 18, 25], the constant state solution in the two-phaseregion 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 evapo-ration rate. A larger evaporation rate leads to faster growth of the disturbances. Oscillationin solutions can occur in the two-phase region before the dry region completely takes over.

We show that the reaction-diffusion 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


Page 4: transportul umiditatii

8/3/2019 transportul umiditatii 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 reaction-diffusion 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 pseudo-spectral method are usedfor spatial discretization.

The rest of the paper is organized as follows. In section 2, we propose a reaction-diffusion

model based on the Hertz-Knudsen 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 reaction-diffusion 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 Reaction-Diffusion 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 







+ λΓ, (1)








− Γ

ρ, (2)










ρ. (3)


Page 5: transportul umiditatii

8/3/2019 transportul umiditatii 5/20

Here Γ is the rate of the phase change (mass per unit volume per unit time), given bythe modified Hertz-Knudsen equation (See [13] where this equation is used for studyingevaporation phenomenon)

Γ = E (1− φ)



(P v − cP s)√T 


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 re-write (4) as

Γ =E (1− φ)ρ



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 pre-heatedat temperature T a by a radiator at temperature T r. The boundary conditions are mixedcondition at x = 0:


∂x= hr(T r − T ) + hc(T a − T ), (7a)


∂x= hv(V a − V ), (7b)



= hw(W a−

W ) (7c)

and symmetric condition at x = L:






∂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.


Page 6: transportul umiditatii

8/3/2019 transportul umiditatii 6/20

2.1 Nondimensionalization

We now proceed to nondimensionalize the equations by choosing the following scaling

θ =T − T 0

T , τ  =


t, ξ =


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


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 param-eters 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)


β  = λE (1− φ)ρL2



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)


Page 7: transportul umiditatii

8/3/2019 transportul umiditatii 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)


h1 = hV L,

h2 = hW L,

h3 =σT 3L


[(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 reaction-diffusion 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 Reaction-diffusion model

Splitting scheme. We adopt a splitting method for the reaction-diffusion model by sepa-rating 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 reaction-diffusion equations and we referinterested readers to [1, 8] and references therein.


Page 8: transportul umiditatii

8/3/2019 transportul umiditatii 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 semi-discrete form can be describedas follows.

1. We solve for vapor content first using

V c

−V (n)

∆τ =−

αS 1


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 non-negative. 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)


Page 9: transportul umiditatii

8/3/2019 transportul umiditatii 9/20

6. Update the temperature using the rate of phase change S 2 given by (23)



∆τ  =


ρ0 + W (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)

∆τ =


ρ0 + W (n+1)θ(n+1)ξξ + S 1 + S 2. (33)


Page 10: transportul umiditatii

8/3/2019 transportul umiditatii 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 reaction-diffusion model (24)-(26) are the same as thatfor the simplified instantaneous model (36)-(38). The only difference between the reaction-diffusion 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 reaction-diffusion 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.


Page 11: transportul umiditatii

8/3/2019 transportul umiditatii 11/20


ν  =∆τ αβ 

√θ + θ0


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 reaction-diffusion 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 reaction-diffusion 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 re-sults 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 non-smoothnessof 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 reaction-diffusion 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 reaction-diffusion 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


Page 12: transportul umiditatii

8/3/2019 transportul umiditatii 12/20

0 20 40 60 80 1000




0 20 40 60 80 1000




0 20 40 60 80 1000




t (min)



0 20 40 60 80 1000




0 20 40 60 80 1000




0 20 40 60 80 1000




t (min)



Figure 1: Computational results of  V , W , and  T  obtained by using the reaction-diffusion 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].


Page 13: transportul umiditatii

8/3/2019 transportul umiditatii 13/20



















 t   (  mi  n )  

x (m)




















 t   (  mi  n )  

x (m)


(a) (b)

Figure 2: Snap shots of the water content computed using the reaction-diffusion 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 dry-wet 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 reaction-diffusion 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ξ).


Page 14: transportul umiditatii

8/3/2019 transportul umiditatii 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


, (42)

sw = −D2m2w + αβ 

−V 1θ + v +

V 0ρ0 + W 0


, (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


y =



, 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.


ρ0 + W 0= θ20D1. (45)

For our choice of parameters in this bread-baking 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


+D2m2(θ20D1m2 + βV 1 + αβ )− αβV 0ρ0 + W 0


= 0. (46)


Page 15: transportul umiditatii

8/3/2019 transportul umiditatii 15/20

3.1.1 Reaction-only

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 Reaction-diffusion

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


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 math-ematical 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 non-positive 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


Page 16: transportul umiditatii

8/3/2019 transportul umiditatii 16/20

0 0.2 0.4 0.6 0.8 1−12









6x 10


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 reaction-diffusion 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.


Page 17: transportul umiditatii

8/3/2019 transportul umiditatii 17/20

0 0.002 0.004 0.006 0.008 0.010




x (cm)

0 0.002 0.004 0.006 0.008 0.010




0 0.002 0.004 0.006 0.008 0.010




0 0.002 0.004 0.006 0.008 0.010




0 0.002 0.004 0.006 0.008 0.010




0 0.002 0.004 0.006 0.008 0.010




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 two-phase region will become unstable. The fastest growing mode eventually causesa multiple region of dry and two-phase 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 ).


Page 18: transportul umiditatii

8/3/2019 transportul umiditatii 18/20

0 20 40 60 80 1000









t (min)

   W  m   i   d ω=0.5




0 2 4 6 8 100









t (min)

   W  m   i   d



(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, pre-vious 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 reaction-diffusion 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 reaction-diffusion model,we showed that the numerical instability can be eliminated by using a sufficiently small timestep. This is due to the fact that the reaction-diffusion 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 two-phase zones.

Our analysis of the reaction-diffusion 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


Page 19: transportul umiditatii

8/3/2019 transportul umiditatii 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 reaction-diffusion 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 R-146-000-053-112 and R-146-000-181-112 from theNational University of Singapore (PL and WZ).


[1] D. Alemani, B. Chopard, J. Galceran, and J. Buffle, LBGK method coupled to timesplitting technique for solving reaction-diffusion processes in complex systems, Phys.Chem. Chem Phys. 7 (2005), no. 18, 3331-3341.

[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), 87-99.

[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, North-Holland,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 reaction-diffusion

systems, Math. Comp. 70 (2000), no. 236, 1481-1501.

[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.174-188.

[10] A.C. Fowler, Mathematical Models in the Applied Sciences, Cambridge University Press,1997.


Page 20: transportul umiditatii

8/3/2019 transportul umiditatii 20/20

[11] M.L. Frankel, G. Kovacic, V. Roytburd and I. timofeyev, Finite-dimensional dynamicalsystem modeling thermal instabilities, Physica D  137 (2000), 295-315.

[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, 1101-1112.

[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), 255-339.

[15] K.A. Landman, C.P. Please, Modelling moisture uptake in a cereal grain, IMA J. Math.Buss. Ind. 10 (2000), 265-287.

[16] M.J. McGuinness, C.P. Please, N. Fowkes, P. McGowan, L. Ryder and D. Forte, Mod-elling the wetting and cooking of a single cereal grain, IMA J. Math. Buss. Ind. 11(2000), no.1, 49-70.

[17] J. D. Murray, Mathematical Biology , Springer-Verlag, 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), 658-663.

[21] K. Thorvaldsson and H. Janestad, A model for simultaneous heat, water and vapordiffusion. J. Food Eng. 40 (1999) 167-172.

[22] Tong, C.H. and Lund, D.B., Microwave heating of baked dough products with simulta-neous heat and moisture transfer, Journal of Food Engineering , 19 (1993), 319-339.

[23] B. Zanoni, C. Peri and S. Pierucci, A study of the bread-baking process I: a phenomeno-

logical model, Journal of Food Engineering , 19 (1993), No. 4, 389-398.

[24] B. Zanoni, C. Peri and S. Pierucci, Study of bread baking process-II: Mathematicalmodelling, Journal of Food Engineering , 23 (1994), No. 3, 321-336.

[25] W Zhou, Application of FDM and FEM to solving the simultaneous heat and mois-ture transfer inside bread during baking, International Journal of Computational Fluid Dynamics, 19 (2005), No. 1, 73-77.