- Email: [email protected]
Mathematical modelling of syneresis of cheese curd
Mathematics and Computers in Simulation 65 (2004) 165–175
Mathematical modelling of syneresis of cheese curd E. Tijskens∗ , J. De Baerdemaeker Laboratory for Agro-Machinery and Agro-Processes, Catholic University Leuven, Kasteelpark Arenberg 30, B-3001 Leuven, Belgium
Abstract A mathematical model for syneresis of cheese curd in one dimension is presented based on the theory of mechanics of porous media. The simple mathematics leave room for surprising physical considerations and several interesting observations are made. Model parameters are estimated and validated using one-dimensional experiments. The agreement between experiment and simulation is reasonable given the speculative nature of the constitutive equations. The model also reproduces the qualitative aspects as observed by MRI experiments. © 2003 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Syneresis; Deforming porous media; Mass transfer
1. Introduction In order to predict and control the process of soft cheese production, cheese-making companies have expressed their interest in a general mathematical model for the drainage of beds of curd grains. From a physical point of view the manufacturing of soft cheeses consists of: (1) coagulating milk to form a gel, the curd, consisting of a protein (casein) skeleton, and an interstitial fluid, the whey, (2) cutting the gel into small pieces (typically 2 cm cubes), (3) pouring the curd grains in a perforated mould to let the whey evacuate from the network. The entire process typically takes about 15 h. The evacuation of the whey is in fact a two step process. The first step is the expulsion of whey by the shrinking curd grains. This process is called syneresis. It is a porous flow with the porous medium defined by the protein skeleton. We refer to this porous medium as the micro-scale porosity. Once the whey is outside the curd grain it enters a porous medium on a larger scale defined by the curd grains themselves. This is referred to as the macro-scale porosity. Here, the channels available for fluid flow are much larger than the channels in the curd grains, but the volume fraction of these channels with respect to the total volume of the curd bed becomes very small soon after the curd grains are poured into the mould. This is because the curd grains are slightly more dense (about 4%) than the whey and they deform readily due to the low strength of ∗
Corresponding author. Tel.: +32-16-32-8595; fax: +32-16-32-8590. E-mail addresses: [email protected] (E. Tijskens), [email protected] (J. De Baerdemaeker). 0378-4754/$30.00 © 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2003.09.016
166
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
the protein skeleton, so they get packed. In addition, at contacts between curd grains there is curd grain fusion [1]. In this paper, we will further concentrate on syneresis, as it is a basic building block in the entire process. The modelling of syneresis of cheese curd is a difficult problem for several reasons. From the physical point of view, the problem is incompletely understood on nearly every scale of interest. On the atomistic scale most work is of physico-chemical nature, and has led to the hairy micelle concept [2,3]. Only rough models of the physics of gel-formation have been developed [4]. The microscopic geometry of the protein skeleton has been revealed by confocal laser scanning images [5] and by transmission electron microscopy [6]. They show a sponge-like structure with channel diameters of about 5 m. A major difficulty of the problem is that this skeleton is not in a stable state. The protein material is continuously rearranging itself towards equilibrium. In doing so, it exerts a pressure on the contained fluid, which responds by trying to escape at the curd grain boundaries while the curd grains continue to shrink, and their porosity decreases. Thus, to some extent, the shrinking curd grain can be regarded as a “self-squeezing sponge”. Most work on syneresis of curd refers to the macro-scale and is experimental [7–11]. Appropriate continuum models for syneresis of curd must treat the curd as a porous medium in which each of the components–that is the protein network and the interstitial fluid—is represented by an overlapping continuum [12–14]. As a porous medium cheese curd is certainly a very unusual material, as it has a very high porosity—initially about 90%—and is very deformable. For example, when put on a flat surface, a curd grain deforms almost instant taneously under its own weight (probably by rupturing of the network). On the other hand, when the curd grain is surrounded by whey, the deformation is slow and almost completely driven by syneresis. However, on the time scale of the entire process, the deformation of curd grains certainly cannot be neglected, as the average porosity evolves from 90 to about 45% implying a volume change of about 50%. Clearly, the problem is not only a problem of fluid flow, but also of deformation of the curd grains. This makes the syneresis problem also a difficult computational problem, as the governing equations will constitute a time dependent nonlinear system of coupled PDEs with a free moving boundary and large deformations, and although the basic equations governing the system at the macro-scale are fairly obvious—conservation of mass and momentum for the whey and the protein network—the selection of constitutive equations describing the material properties is subject to much uncertainty. Now that the general system and its difficulties are sketched, the paper will proceed as follows. In view of all the uncertainties on constitutive equations associated with the problem, it was decided that a model for syneresis should first consider a simplified situation before the problem should be attacked in three-dimensional. Thus, a one-dimensional syneresis problem is sketched Section 2 and a basic mathematical model is developed. The model is simple, yet it identifies accurately the physical properties that come into play and that are needed in a three-dimensional model. Its development has increased our insight in the syneresis problem significantly. Section 3 describes an experimental realisation of this model. It is used to estimate the material properties, and to validate the model. The experiments show that the basic model requires some adaptation and a physical explanation for this is provided. Finally, Section 4 summarises the conclusions.
2. Syneresis modelling in one-dimensional Consider the following problem in which thin cylindrical layers of curd are prepared in a large petri disk (Fig. 1). Because the curd is prepared in the petri disk it is attached to the walls, the whey can only
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
167
z 1D behaviour h (t=0) h (t>0) 0
Fig. 1. A one-dimensional syneresis problem.
be expelled at the top surface. In the centre of the disk the fluid flow is directed along the vertical axis because of symmetry. Consequently, the syneresis is one-dimensional at the centre of the disk. According to Coussy’s treatment of the mechanics of porous media [12], the curd grain is treated as a system of two overlapping continua, the skeleton, representing the protein network, and the fluid, representing the whey. Here, both components are considered here as ideal incompressible fluids, i.e. they are inviscid. Hence, the mass balance for the fluid, respectively the skeleton reads: ∂φ ∂φV ∂qfl + + = 0, ∂t ∂z ∂z −
∂φ ∂ (1 − φ) V + = 0. ∂t ∂z
(1) (2)
Here, φ is the local (micro-scale) porosity, V the velocity of the skeleton particle, qfl the volume flux of the fluid relative to the skeleton, t represents time and z is a spatial co-ordinate along the vertical axis, pointing upwards. In Eqs. (1) and (2), the constant densities have been divided away. Summing both expressions yields the total mass balance: ∂V ∂qfl + = 0. ∂z ∂z
(3)
This equation expresses that the skeleton has to move at the same speed as the fluid, but in the opposite direction. This is a direct consequence of the assumption that both components are incompressible. The slow movement of fluid allows for modelling the fluid flow by Darcy’s law: ∂p fl fl q =K − −ρ g . (4) ∂z In Eq. (4), the permeability K is modelled by the Carman–Kozeny equation: K(φ) = κ
φ3 d2 . (1 − φ)2 m
(5)
Here, dm is the mean channel diameter, and κ depends on the microstructure. For constant κ and dm , the permeability K(φ) is monotonously increasing with the porosity φ. As the curd grain shrinks, the porosity obviously decreases, and thinking of the curd grain as a “self-squeezing sponge”, one would also expect the mean channel diameter dm to decrease, so that the permeability also decreases. This is not the case, however, as is shown in experiments with curd by Dejmek [15]. It was found that constrained curd—that is curd which is practically not allowed to shrink because it is attached to the wall—shows an increasing permeability with time. As the porosity of constrained curd is almost constant, this suggests that dm must
168
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
be increasing with time. This is another consequence of the dynamical rearrangement of the skeleton towards an equilibrium state. Thus, we assume: dm (t) = D(t)dm (0), D(t) = 1 + D∞ (1 − e−αt ),
(6)
where D(t) satisfies D(t = 0) = 1. If we absorb the constant factor dm (0)2 in k we obtain K(φ) = k
φ3 D(t)2 . (1 − φ)2
(7)
We proceed by decomposing the fluid pressure inside the curd grain in a hydrostatic component ph and a syneresis component or overpressure ps : p = ps + ph , ph = p(z1 ) −
z z1
gρ(z) dz ,
where the reference level z1 is taken to be the curd–whey interface. Hence, we obtain ∂ps fl fl + (ρ − ρ )g . q =K − ∂z
(8)
(9)
It is unclear which density one has to take for ρ. Two interpretations are possible. Either the fluid does not feel the weight of the skeleton and its hydrostatic pressure is caused by its own weight only, in which case ρ = ρfl and the second term in Eq. (9) vanishes: qfl = −K
∂ps . ∂z
At the other extreme, the fluid does feel the weight of the skeleton, i.e. ρ = ρcurd , and ∂ps fl + (1 − φ)ρg . q =K − ∂z
(10)
(11)
Here, we defined ρ ≡ ρsk − ρfl ≈ 0.1 kg/m3 .
(12)
The extra driving force fwsk = (1 − φ)ρg
(13)
is identical to the driving force determining the settling of a system of fluid particles in suspension. We shall therefor refer to it as the “settling” force. The obvious time evolution of a settling system is that the heavy particles accumulate at the bottom while the lighter interstitial fluid is pushed upwards. The first case, described by Eq. (10), can be recovered from Eq. (11) by setting ∆ρ = 0. An argument in favour of the case ρ = ρfl , is that a curd grain submerged in whey can be observed not to collapse under its own weight, as it does in free air. The observed deformation of a curd grain submerged in whey is due to syneresis (shrinking). The computations (see below) show that the magnitude of this driving force generally exceeds the syneresis component by two orders of magnitude, except where
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
169
the syneresis pressure gradient, ∂ps /∂z, is very steep. The inclusion of the gravity term in simulations accordingly leads to unrealistically high predicted syneresis rates. It may, however, be argued that this is due to the skeleton being treated as an inviscid fluid, while in reality its viscosity exceeds that of the fluid by several orders of magnitude, causing a finite resistance to deformation and balancing the gravity forces. That is, the settling force is compensated by the viscous friction of the skeleton and either both or none have to be included. For the syneresis pressure ps , the following assumption seems acceptable: ps = β
δv − δv∞ δv
(14)
for an infinitesimally small volume δv. Here, δv∞ is the equilibrium value of δv. By taking account of the definition of porosity and the incompressibility assumption Eq. (14) can be expressed in terms of porosity: ps = β
φ − φ∞ 1−φ
(15)
Finally, integrating Eq. (3) to obtain V(z), and substituting the result into Eq. (1) we obtain a scalar partial differential equation in three-dimensional: ∂φ ∂φ ∂(1 − φ)qfl + qfl (z0 ) + =0 ∂t ∂z ∂z
(16)
where z0 is an arbitrary reference level. It should be noted that Eq. (3) can only be integrated in the one-dimensional case. In two or more dimensions, Eq. (3) can no longer be integrated to obtain V(z), and a momentum balance becomes necessary to close the system of equations. Eq. (16) is subject to boundary conditions φ(z1 , t > 0) = φ∞
(17)
qfl (z = 0, t > 0) = 0
(18)
and initial condition φ(0 ≤ z ≤ z1 , t = 0) = φ0
(19)
By letting the reference level z0 coincide with z = 0, which is the impermeable boundary (qfl = 0) at the bottom of the petri disk, the middle term in Eq. (16) vanishes. The partial differential Eqs. (16)–(19) is discretised using a standard Galerkin finite element method (see, e.g. [16]). As the problem is nonlinear it is transformed into a series of linear PDEs by means of a Newton–Raphson scheme [17]. The Jacobians needed for the Newton–Raphson scheme are computed exactly and efficiently by automatic differentiation using Fastder++ [18,19]. The typical time evolution of the porosity profile along the central axis of the experimental set-up (see Fig. 1) obtained with this model is shown in Fig. 2, without inclusion of the settling force, and in Fig. 3, with inclusion of the settling force. The difference between both is obvious. As the settling force dominates the syneresis pressure by approximately two orders of magnitude, settling is the dominant process, yielding a positive slope throughout most of the domain, and an almost vertical slope at the curd–whey interface (Fig. 3). Thus, the porosity decreases towards the bottom. Without settling force, the
170
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
Fig. 2. Typical time evolution of the porosity profile when the settling force is neglected (problem parameters: h(t = 0) = 0.02, kβ = 10−10 , φ0 = 0.9, φ∞ = 0.5, H(z) = D(t) = 1, all quantities in SI units).
Fig. 3. Typical time evolution of porosity profile with the settling force included (same problem parameters as in Fig. 2).
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
171
Fig. 4. MRI image of a system of syneresing curd grains after 2 hours and 37 minutes. The bottom of the mould is perforated. A recipient is positioned below the mould. Curd grain contours are visible as low humidity areas. The arrows show the direction of fluid transport at the boundary.
profile shows a decrease towards the top of the curd–whey interface and a much slower time evolution. Only the latter is in agreement with MRI observations [20]. Fig. 4 clearly shows that syneresing curd grains develop a decreasing porosity towards the curd–whey interface. Curd grain contours are visible as low humidity areas. These low humidity areas seem to develop rather uniformly all around the grains showing that the porosity profile across a curd grain is quite independent of the orientation of the interface with respect to the vertical axis. 3. Parameter estimation and validation Experiments of a one-dimensional model system for syneresis of curd have been carried out by Lodaite et al. [21]. Thin layers of curd (initial height, h0 = 3, 6.5, 10.5, 20 and 35 mm) were prepared in a large petri disk. After preparation the curd surface is dry, and surface tension counteracts the syneresis pressure. As soon as the surface is wetted, syneresis proceeds. When the surface is not wetted, the skeleton will still rearrange itself, but at a constant porosity. This process is referred to as ageing. The height of the curd layer was measured by means of a laser displacement sensor (see Fig. 1). The shrinkage, defined as h(t) − h(t = 0)
(20)
is plotted in Fig. 5. In one dimension, the shrinkage is related to the porosity φ by the volume balance: h(t) h(0) h(t) − h(0) = φ(t) dh − φ(0) dh . (21) 0
0
Eq. (16) is entirely local in nature and therefore predicts the initial behaviour to be independent of the thickness of the curd layer. This is demonstrated in Fig. 6, where the slope of all curves at time t = 0
172
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
shrinkage [m]
0
x 10
-3
-1
0.0030
-2
0.0065
-3
0.0105
-4
0.0200 -5
0.0350 -6 0
1000
2000
3000
4000 time [s]
5000
6000
7000
8000
Fig. 5. Experimental shrinkage as defined by Eq. (20) of the curd layer as a function of time for different initial heights (m), h(t = 0) = {0.003, 0.0065, 0.0105, 0.020 and 0.035} at pH = 6.0.
is identical. The experiments, however, shown in Fig. 5, show a quite different picture. There, the slope of the curve at time t = 0 clearly varies systematically with the initial height of the curd bed h(t = 0). This thickness dependence is thought to originate from the fact that the skeleton is attached to the bottom of the petri disk and therefore cannot pursue its normal time evolution. This is called the wall effect. To -3
0
h040
x 10
-1 -2
shrinkage h-h0 []
-3 -4 -5 -6 -7 -8
0
1000
2000
3000
4000 time [s]
5000
6000
7000
8000
Fig. 6. Effect of the initial height on the shrinkage h(t) − h0 in the absence of a wall effect (H(z) = 1). From top to bottom: h(t = 0) = {0.01, 0.02, 0.04 and 0.08} (kβ = 5 × 10−11 , other model parameters as in Fig. 2).
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
173
Table 1 Estimated model parameters for the model described by Eqs. (16)–(19) and (22), accounting for the wall effect (all quantities in SI units) kβ = 2.20 × 10−11 D∞ = 0.304 H0 = 0.0989
φ∞ = −0.00699 α = 4.38 × 10−4 λ = 88.6
account for such a wall effect, the constitutive equation for the permeability was adapted as φ3 D(t)2 (1 − φ)2 H(z) = H0 + (1 − H0 )(1 − e−λz )
K = H(z)κ
(22)
The effect of the function H(z) is to diminish the increase of the mean channel diameter towards the bottom of the petri disk, where the skeleton is attached to the glass wall. The model parameters (see Table 1) have been estimated indirectly for two of the five experiments (h(t = 0) = 0.00105 and 0.0020) by minimising the error function 2 (hsim,h0 =0.0105 (ti ) − hexp,h0 =0.0105 (ti )) + (23) (hsim,h0 =0.020 (ti ) − hexp,h0 =0.020 (ti ))2 i
i
The remaining three experiments (h(t = 0) = 0.0030, 0.0065 and 0.035) are used for validation. The results are shown in Fig. 7. They are remarkably good considering the speculative nature of the constitutive
shrinkageh-h0[m]
0
x 10
-3
-1
0.0030
-2
0.0065
-3
0.0105
-4
0.0200 -5
0.0350 -6
0
1000
2000
3000
4000 time [s]
5000
6000
7000
8000
Fig. 7. Comparison of experimental (see Fig. 5) and simulated shrinkage curves. The simulated curves (dashed lines) have been obtained using the parameters of Table 1. The curves with initial heights h(t = 0) = {0.00105 and 0.0020} have been used for parameter estimation. The remaining curves for validation.
174
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
equations and the simplicity of the model. It also shows that the difference in syneresis rate for different initial thicknesses is probably due to the external influence of the bottom wall rather than due to a nonlocal driving force. The slightly negative value for φ∞ is considered an indication that φ∞ should in fact equal zero, although in practice the porosity is never observed to approach zero. The reason for this is that, as time goes by, the skeleton becomes firmer due to ageing, and is no longer able to rearrange. The function H(z) approaches unity as z goes to infinity. In practice, using the estimated parameters from Table 1, H(z) ∼ 0.95 for z = 0.035 m. Thus, the wall effect is estimated to penetrate 3.5–4.5 cm into the curd. 4. Conclusions The one-dimensional syneresis of curd can be modelled by considering the skeleton and the fluid as an ideal fluid (incompressible and inviscid). The fact that the skeleton actually has a very high viscosity, is compensated by neglecting the settling force fwsk (13). The governing Eqs. (16)–(19) permit a reasonably good prediction of the experiments, given the speculative nature of the constitutive equations for the permeability (7) and (22) and the syneresis pressure (15). If the curd is attached to the wall, this has a significant effect on the ability of the skeleton to dynamically rearrange and on the syneresis process. The penetration depth of this wall effect is estimated to be in the range 0.035–0.045 m. Despite the simplicity of the model, it clearly identifies all the important properties that will be needed to develop a full three-dimensional model. The model is in fact overly simple, yet it accurately identifies the physical properties that come into play and that will be needed in a full three-dimensional model of syneresing curd grains. A first improvement would be to treat the skeleton as a viscous fluid, and incorporate the settling force. Acknowledgements This research was sponsored by the European Community under contract Fair-CT-1096. The authors acknowledge Professor Peter Dejmek and Dr. Karin Östergren of Lund University (Sweden) and Dr. François Mariette of Cemagref Rennes (France) for numerous interesting discussions and for providing the experimental data. References [1] [2] [3] [4] [5] [6] [7] [8]
J.C. Akkerman, R.O. Lewis, P. Walstra, Fusion of curd grains, Neth. Milk Dairy J. 47 (1993) 137–144. P. Walstra, T. van Vliet, The physical chemistry of curd making, Neth. Milk Dairy J. 40 (1986) 241–259. C. Holt, D.S. Horne, The hairy casein micelle: evolution of the concept, Neth. Milk Dairy J. 50 (1996) 85–111. M.T.A. Bos, J.H.J. Opheusden, Brownian dynamics simulation of gelation and ageing in interacting systems, Phys. Rev. E 53 (5) (1996) 5044–5050. A.N. Hassan, J.F. Frank, Modification of microstructure and texture of rennet curd by using a capsule-forming non-ropy lactic culture, J. Dairy Res. 64 (1997) 115–121. J.M. Aguilera, J.E. Kinsella, Compression strength of dairy gels and microstructural interpretation, J. Food Sci. 56 (5) (1991) 1224–1228. H.J.M. van Dijk, Syneresis of curd, Ph.D. thesis, Landbouwuniversiteit Wageningen, Wageningen, The Netherlands, 1982. P. Zoon, Rheological properties of rennet-induced skim milk gels, Ph.D. thesis, Landbouwuniversiteit Wageningen, Wageningen, The Netherlands, 1988.
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
175
[9] J.C. Akkerman, Drainage of curd, Ph.D. thesis, Landbouwuniversiteit Wageningen, Wageningen, The Netherlands, 1992. [10] J.C. Akkerman, F.H.J. Fox, P. Walstra, Drainage of curd: expression of single curd grains, Neth. Milk Dairy J. 48 (1994) 1–17. [11] J.C. Akkerman, F.H.J. Fox, P. Walstra, Drainage of curd: role of drainage equipment in relation to curd properties, Neth. Milk Dairy J. 50 (1996) 371–406. [12] O. Coussy, Mechanics of Porous Continua, Wiley, New York, 1995, p. 454. [13] R.W. Lewis, B.A. Schrefler, The Finite Element Method in the Deformation and Consolidation of Porous Media. Wiley, New York, 1998. [14] R. DeBoer, Theory of Porous Media: Highlights in Historical Development and Current State, Springer-Verlag, Berlin, 1999. [15] P. Dejmek, Personal communication, 1999. [16] D.S. Burnett, Finite Element Analysis, Addison-Wesley, Reading, MA, 1987. [17] C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, PA, 1995. [18] E. Tijskens, D. Roose, H. Ramon, J. De Baerdemaeker, Automatic differentiation for solving nonlinear partial differential equations: an efficient operator overloading approach, Numer. Algorith. 30 (2002) 259–301. [19] E. Tijskens, H. Ramon, J. De Baerdemaeker, FastDer++: efficient automatic differentiation for nonlinear PDE solvers, Math. Comput. Simul. 65 (2004) 177–190. [20] F. Mariette, Personal communication, 1998. [21] C. Lodaite, K. Östergren, M. Paulsson, P. Dejmek, One-dimensional syneresis of rennet-induced gels, Int. Dairy J. 10 (12) (2000) 829–834.
Mathematical modelling of syneresis of cheese curd E. Tijskens∗ , J. De Baerdemaeker Laboratory for Agro-Machinery and Agro-Processes, Catholic University Leuven, Kasteelpark Arenberg 30, B-3001 Leuven, Belgium
Abstract A mathematical model for syneresis of cheese curd in one dimension is presented based on the theory of mechanics of porous media. The simple mathematics leave room for surprising physical considerations and several interesting observations are made. Model parameters are estimated and validated using one-dimensional experiments. The agreement between experiment and simulation is reasonable given the speculative nature of the constitutive equations. The model also reproduces the qualitative aspects as observed by MRI experiments. © 2003 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Syneresis; Deforming porous media; Mass transfer
1. Introduction In order to predict and control the process of soft cheese production, cheese-making companies have expressed their interest in a general mathematical model for the drainage of beds of curd grains. From a physical point of view the manufacturing of soft cheeses consists of: (1) coagulating milk to form a gel, the curd, consisting of a protein (casein) skeleton, and an interstitial fluid, the whey, (2) cutting the gel into small pieces (typically 2 cm cubes), (3) pouring the curd grains in a perforated mould to let the whey evacuate from the network. The entire process typically takes about 15 h. The evacuation of the whey is in fact a two step process. The first step is the expulsion of whey by the shrinking curd grains. This process is called syneresis. It is a porous flow with the porous medium defined by the protein skeleton. We refer to this porous medium as the micro-scale porosity. Once the whey is outside the curd grain it enters a porous medium on a larger scale defined by the curd grains themselves. This is referred to as the macro-scale porosity. Here, the channels available for fluid flow are much larger than the channels in the curd grains, but the volume fraction of these channels with respect to the total volume of the curd bed becomes very small soon after the curd grains are poured into the mould. This is because the curd grains are slightly more dense (about 4%) than the whey and they deform readily due to the low strength of ∗
Corresponding author. Tel.: +32-16-32-8595; fax: +32-16-32-8590. E-mail addresses: [email protected] (E. Tijskens), [email protected] (J. De Baerdemaeker). 0378-4754/$30.00 © 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2003.09.016
166
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
the protein skeleton, so they get packed. In addition, at contacts between curd grains there is curd grain fusion [1]. In this paper, we will further concentrate on syneresis, as it is a basic building block in the entire process. The modelling of syneresis of cheese curd is a difficult problem for several reasons. From the physical point of view, the problem is incompletely understood on nearly every scale of interest. On the atomistic scale most work is of physico-chemical nature, and has led to the hairy micelle concept [2,3]. Only rough models of the physics of gel-formation have been developed [4]. The microscopic geometry of the protein skeleton has been revealed by confocal laser scanning images [5] and by transmission electron microscopy [6]. They show a sponge-like structure with channel diameters of about 5 m. A major difficulty of the problem is that this skeleton is not in a stable state. The protein material is continuously rearranging itself towards equilibrium. In doing so, it exerts a pressure on the contained fluid, which responds by trying to escape at the curd grain boundaries while the curd grains continue to shrink, and their porosity decreases. Thus, to some extent, the shrinking curd grain can be regarded as a “self-squeezing sponge”. Most work on syneresis of curd refers to the macro-scale and is experimental [7–11]. Appropriate continuum models for syneresis of curd must treat the curd as a porous medium in which each of the components–that is the protein network and the interstitial fluid—is represented by an overlapping continuum [12–14]. As a porous medium cheese curd is certainly a very unusual material, as it has a very high porosity—initially about 90%—and is very deformable. For example, when put on a flat surface, a curd grain deforms almost instant taneously under its own weight (probably by rupturing of the network). On the other hand, when the curd grain is surrounded by whey, the deformation is slow and almost completely driven by syneresis. However, on the time scale of the entire process, the deformation of curd grains certainly cannot be neglected, as the average porosity evolves from 90 to about 45% implying a volume change of about 50%. Clearly, the problem is not only a problem of fluid flow, but also of deformation of the curd grains. This makes the syneresis problem also a difficult computational problem, as the governing equations will constitute a time dependent nonlinear system of coupled PDEs with a free moving boundary and large deformations, and although the basic equations governing the system at the macro-scale are fairly obvious—conservation of mass and momentum for the whey and the protein network—the selection of constitutive equations describing the material properties is subject to much uncertainty. Now that the general system and its difficulties are sketched, the paper will proceed as follows. In view of all the uncertainties on constitutive equations associated with the problem, it was decided that a model for syneresis should first consider a simplified situation before the problem should be attacked in three-dimensional. Thus, a one-dimensional syneresis problem is sketched Section 2 and a basic mathematical model is developed. The model is simple, yet it identifies accurately the physical properties that come into play and that are needed in a three-dimensional model. Its development has increased our insight in the syneresis problem significantly. Section 3 describes an experimental realisation of this model. It is used to estimate the material properties, and to validate the model. The experiments show that the basic model requires some adaptation and a physical explanation for this is provided. Finally, Section 4 summarises the conclusions.
2. Syneresis modelling in one-dimensional Consider the following problem in which thin cylindrical layers of curd are prepared in a large petri disk (Fig. 1). Because the curd is prepared in the petri disk it is attached to the walls, the whey can only
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
167
z 1D behaviour h (t=0) h (t>0) 0
Fig. 1. A one-dimensional syneresis problem.
be expelled at the top surface. In the centre of the disk the fluid flow is directed along the vertical axis because of symmetry. Consequently, the syneresis is one-dimensional at the centre of the disk. According to Coussy’s treatment of the mechanics of porous media [12], the curd grain is treated as a system of two overlapping continua, the skeleton, representing the protein network, and the fluid, representing the whey. Here, both components are considered here as ideal incompressible fluids, i.e. they are inviscid. Hence, the mass balance for the fluid, respectively the skeleton reads: ∂φ ∂φV ∂qfl + + = 0, ∂t ∂z ∂z −
∂φ ∂ (1 − φ) V + = 0. ∂t ∂z
(1) (2)
Here, φ is the local (micro-scale) porosity, V the velocity of the skeleton particle, qfl the volume flux of the fluid relative to the skeleton, t represents time and z is a spatial co-ordinate along the vertical axis, pointing upwards. In Eqs. (1) and (2), the constant densities have been divided away. Summing both expressions yields the total mass balance: ∂V ∂qfl + = 0. ∂z ∂z
(3)
This equation expresses that the skeleton has to move at the same speed as the fluid, but in the opposite direction. This is a direct consequence of the assumption that both components are incompressible. The slow movement of fluid allows for modelling the fluid flow by Darcy’s law: ∂p fl fl q =K − −ρ g . (4) ∂z In Eq. (4), the permeability K is modelled by the Carman–Kozeny equation: K(φ) = κ
φ3 d2 . (1 − φ)2 m
(5)
Here, dm is the mean channel diameter, and κ depends on the microstructure. For constant κ and dm , the permeability K(φ) is monotonously increasing with the porosity φ. As the curd grain shrinks, the porosity obviously decreases, and thinking of the curd grain as a “self-squeezing sponge”, one would also expect the mean channel diameter dm to decrease, so that the permeability also decreases. This is not the case, however, as is shown in experiments with curd by Dejmek [15]. It was found that constrained curd—that is curd which is practically not allowed to shrink because it is attached to the wall—shows an increasing permeability with time. As the porosity of constrained curd is almost constant, this suggests that dm must
168
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
be increasing with time. This is another consequence of the dynamical rearrangement of the skeleton towards an equilibrium state. Thus, we assume: dm (t) = D(t)dm (0), D(t) = 1 + D∞ (1 − e−αt ),
(6)
where D(t) satisfies D(t = 0) = 1. If we absorb the constant factor dm (0)2 in k we obtain K(φ) = k
φ3 D(t)2 . (1 − φ)2
(7)
We proceed by decomposing the fluid pressure inside the curd grain in a hydrostatic component ph and a syneresis component or overpressure ps : p = ps + ph , ph = p(z1 ) −
z z1
gρ(z) dz ,
where the reference level z1 is taken to be the curd–whey interface. Hence, we obtain ∂ps fl fl + (ρ − ρ )g . q =K − ∂z
(8)
(9)
It is unclear which density one has to take for ρ. Two interpretations are possible. Either the fluid does not feel the weight of the skeleton and its hydrostatic pressure is caused by its own weight only, in which case ρ = ρfl and the second term in Eq. (9) vanishes: qfl = −K
∂ps . ∂z
At the other extreme, the fluid does feel the weight of the skeleton, i.e. ρ = ρcurd , and ∂ps fl + (1 − φ)ρg . q =K − ∂z
(10)
(11)
Here, we defined ρ ≡ ρsk − ρfl ≈ 0.1 kg/m3 .
(12)
The extra driving force fwsk = (1 − φ)ρg
(13)
is identical to the driving force determining the settling of a system of fluid particles in suspension. We shall therefor refer to it as the “settling” force. The obvious time evolution of a settling system is that the heavy particles accumulate at the bottom while the lighter interstitial fluid is pushed upwards. The first case, described by Eq. (10), can be recovered from Eq. (11) by setting ∆ρ = 0. An argument in favour of the case ρ = ρfl , is that a curd grain submerged in whey can be observed not to collapse under its own weight, as it does in free air. The observed deformation of a curd grain submerged in whey is due to syneresis (shrinking). The computations (see below) show that the magnitude of this driving force generally exceeds the syneresis component by two orders of magnitude, except where
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
169
the syneresis pressure gradient, ∂ps /∂z, is very steep. The inclusion of the gravity term in simulations accordingly leads to unrealistically high predicted syneresis rates. It may, however, be argued that this is due to the skeleton being treated as an inviscid fluid, while in reality its viscosity exceeds that of the fluid by several orders of magnitude, causing a finite resistance to deformation and balancing the gravity forces. That is, the settling force is compensated by the viscous friction of the skeleton and either both or none have to be included. For the syneresis pressure ps , the following assumption seems acceptable: ps = β
δv − δv∞ δv
(14)
for an infinitesimally small volume δv. Here, δv∞ is the equilibrium value of δv. By taking account of the definition of porosity and the incompressibility assumption Eq. (14) can be expressed in terms of porosity: ps = β
φ − φ∞ 1−φ
(15)
Finally, integrating Eq. (3) to obtain V(z), and substituting the result into Eq. (1) we obtain a scalar partial differential equation in three-dimensional: ∂φ ∂φ ∂(1 − φ)qfl + qfl (z0 ) + =0 ∂t ∂z ∂z
(16)
where z0 is an arbitrary reference level. It should be noted that Eq. (3) can only be integrated in the one-dimensional case. In two or more dimensions, Eq. (3) can no longer be integrated to obtain V(z), and a momentum balance becomes necessary to close the system of equations. Eq. (16) is subject to boundary conditions φ(z1 , t > 0) = φ∞
(17)
qfl (z = 0, t > 0) = 0
(18)
and initial condition φ(0 ≤ z ≤ z1 , t = 0) = φ0
(19)
By letting the reference level z0 coincide with z = 0, which is the impermeable boundary (qfl = 0) at the bottom of the petri disk, the middle term in Eq. (16) vanishes. The partial differential Eqs. (16)–(19) is discretised using a standard Galerkin finite element method (see, e.g. [16]). As the problem is nonlinear it is transformed into a series of linear PDEs by means of a Newton–Raphson scheme [17]. The Jacobians needed for the Newton–Raphson scheme are computed exactly and efficiently by automatic differentiation using Fastder++ [18,19]. The typical time evolution of the porosity profile along the central axis of the experimental set-up (see Fig. 1) obtained with this model is shown in Fig. 2, without inclusion of the settling force, and in Fig. 3, with inclusion of the settling force. The difference between both is obvious. As the settling force dominates the syneresis pressure by approximately two orders of magnitude, settling is the dominant process, yielding a positive slope throughout most of the domain, and an almost vertical slope at the curd–whey interface (Fig. 3). Thus, the porosity decreases towards the bottom. Without settling force, the
170
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
Fig. 2. Typical time evolution of the porosity profile when the settling force is neglected (problem parameters: h(t = 0) = 0.02, kβ = 10−10 , φ0 = 0.9, φ∞ = 0.5, H(z) = D(t) = 1, all quantities in SI units).
Fig. 3. Typical time evolution of porosity profile with the settling force included (same problem parameters as in Fig. 2).
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
171
Fig. 4. MRI image of a system of syneresing curd grains after 2 hours and 37 minutes. The bottom of the mould is perforated. A recipient is positioned below the mould. Curd grain contours are visible as low humidity areas. The arrows show the direction of fluid transport at the boundary.
profile shows a decrease towards the top of the curd–whey interface and a much slower time evolution. Only the latter is in agreement with MRI observations [20]. Fig. 4 clearly shows that syneresing curd grains develop a decreasing porosity towards the curd–whey interface. Curd grain contours are visible as low humidity areas. These low humidity areas seem to develop rather uniformly all around the grains showing that the porosity profile across a curd grain is quite independent of the orientation of the interface with respect to the vertical axis. 3. Parameter estimation and validation Experiments of a one-dimensional model system for syneresis of curd have been carried out by Lodaite et al. [21]. Thin layers of curd (initial height, h0 = 3, 6.5, 10.5, 20 and 35 mm) were prepared in a large petri disk. After preparation the curd surface is dry, and surface tension counteracts the syneresis pressure. As soon as the surface is wetted, syneresis proceeds. When the surface is not wetted, the skeleton will still rearrange itself, but at a constant porosity. This process is referred to as ageing. The height of the curd layer was measured by means of a laser displacement sensor (see Fig. 1). The shrinkage, defined as h(t) − h(t = 0)
(20)
is plotted in Fig. 5. In one dimension, the shrinkage is related to the porosity φ by the volume balance: h(t) h(0) h(t) − h(0) = φ(t) dh − φ(0) dh . (21) 0
0
Eq. (16) is entirely local in nature and therefore predicts the initial behaviour to be independent of the thickness of the curd layer. This is demonstrated in Fig. 6, where the slope of all curves at time t = 0
172
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
shrinkage [m]
0
x 10
-3
-1
0.0030
-2
0.0065
-3
0.0105
-4
0.0200 -5
0.0350 -6 0
1000
2000
3000
4000 time [s]
5000
6000
7000
8000
Fig. 5. Experimental shrinkage as defined by Eq. (20) of the curd layer as a function of time for different initial heights (m), h(t = 0) = {0.003, 0.0065, 0.0105, 0.020 and 0.035} at pH = 6.0.
is identical. The experiments, however, shown in Fig. 5, show a quite different picture. There, the slope of the curve at time t = 0 clearly varies systematically with the initial height of the curd bed h(t = 0). This thickness dependence is thought to originate from the fact that the skeleton is attached to the bottom of the petri disk and therefore cannot pursue its normal time evolution. This is called the wall effect. To -3
0
h040
x 10
-1 -2
shrinkage h-h0 []
-3 -4 -5 -6 -7 -8
0
1000
2000
3000
4000 time [s]
5000
6000
7000
8000
Fig. 6. Effect of the initial height on the shrinkage h(t) − h0 in the absence of a wall effect (H(z) = 1). From top to bottom: h(t = 0) = {0.01, 0.02, 0.04 and 0.08} (kβ = 5 × 10−11 , other model parameters as in Fig. 2).
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
173
Table 1 Estimated model parameters for the model described by Eqs. (16)–(19) and (22), accounting for the wall effect (all quantities in SI units) kβ = 2.20 × 10−11 D∞ = 0.304 H0 = 0.0989
φ∞ = −0.00699 α = 4.38 × 10−4 λ = 88.6
account for such a wall effect, the constitutive equation for the permeability was adapted as φ3 D(t)2 (1 − φ)2 H(z) = H0 + (1 − H0 )(1 − e−λz )
K = H(z)κ
(22)
The effect of the function H(z) is to diminish the increase of the mean channel diameter towards the bottom of the petri disk, where the skeleton is attached to the glass wall. The model parameters (see Table 1) have been estimated indirectly for two of the five experiments (h(t = 0) = 0.00105 and 0.0020) by minimising the error function 2 (hsim,h0 =0.0105 (ti ) − hexp,h0 =0.0105 (ti )) + (23) (hsim,h0 =0.020 (ti ) − hexp,h0 =0.020 (ti ))2 i
i
The remaining three experiments (h(t = 0) = 0.0030, 0.0065 and 0.035) are used for validation. The results are shown in Fig. 7. They are remarkably good considering the speculative nature of the constitutive
shrinkageh-h0[m]
0
x 10
-3
-1
0.0030
-2
0.0065
-3
0.0105
-4
0.0200 -5
0.0350 -6
0
1000
2000
3000
4000 time [s]
5000
6000
7000
8000
Fig. 7. Comparison of experimental (see Fig. 5) and simulated shrinkage curves. The simulated curves (dashed lines) have been obtained using the parameters of Table 1. The curves with initial heights h(t = 0) = {0.00105 and 0.0020} have been used for parameter estimation. The remaining curves for validation.
174
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
equations and the simplicity of the model. It also shows that the difference in syneresis rate for different initial thicknesses is probably due to the external influence of the bottom wall rather than due to a nonlocal driving force. The slightly negative value for φ∞ is considered an indication that φ∞ should in fact equal zero, although in practice the porosity is never observed to approach zero. The reason for this is that, as time goes by, the skeleton becomes firmer due to ageing, and is no longer able to rearrange. The function H(z) approaches unity as z goes to infinity. In practice, using the estimated parameters from Table 1, H(z) ∼ 0.95 for z = 0.035 m. Thus, the wall effect is estimated to penetrate 3.5–4.5 cm into the curd. 4. Conclusions The one-dimensional syneresis of curd can be modelled by considering the skeleton and the fluid as an ideal fluid (incompressible and inviscid). The fact that the skeleton actually has a very high viscosity, is compensated by neglecting the settling force fwsk (13). The governing Eqs. (16)–(19) permit a reasonably good prediction of the experiments, given the speculative nature of the constitutive equations for the permeability (7) and (22) and the syneresis pressure (15). If the curd is attached to the wall, this has a significant effect on the ability of the skeleton to dynamically rearrange and on the syneresis process. The penetration depth of this wall effect is estimated to be in the range 0.035–0.045 m. Despite the simplicity of the model, it clearly identifies all the important properties that will be needed to develop a full three-dimensional model. The model is in fact overly simple, yet it accurately identifies the physical properties that come into play and that will be needed in a full three-dimensional model of syneresing curd grains. A first improvement would be to treat the skeleton as a viscous fluid, and incorporate the settling force. Acknowledgements This research was sponsored by the European Community under contract Fair-CT-1096. The authors acknowledge Professor Peter Dejmek and Dr. Karin Östergren of Lund University (Sweden) and Dr. François Mariette of Cemagref Rennes (France) for numerous interesting discussions and for providing the experimental data. References [1] [2] [3] [4] [5] [6] [7] [8]
J.C. Akkerman, R.O. Lewis, P. Walstra, Fusion of curd grains, Neth. Milk Dairy J. 47 (1993) 137–144. P. Walstra, T. van Vliet, The physical chemistry of curd making, Neth. Milk Dairy J. 40 (1986) 241–259. C. Holt, D.S. Horne, The hairy casein micelle: evolution of the concept, Neth. Milk Dairy J. 50 (1996) 85–111. M.T.A. Bos, J.H.J. Opheusden, Brownian dynamics simulation of gelation and ageing in interacting systems, Phys. Rev. E 53 (5) (1996) 5044–5050. A.N. Hassan, J.F. Frank, Modification of microstructure and texture of rennet curd by using a capsule-forming non-ropy lactic culture, J. Dairy Res. 64 (1997) 115–121. J.M. Aguilera, J.E. Kinsella, Compression strength of dairy gels and microstructural interpretation, J. Food Sci. 56 (5) (1991) 1224–1228. H.J.M. van Dijk, Syneresis of curd, Ph.D. thesis, Landbouwuniversiteit Wageningen, Wageningen, The Netherlands, 1982. P. Zoon, Rheological properties of rennet-induced skim milk gels, Ph.D. thesis, Landbouwuniversiteit Wageningen, Wageningen, The Netherlands, 1988.
E. Tijskens, J. De Baerdemaeker / Mathematics and Computers in Simulation 65 (2004) 165–175
175
[9] J.C. Akkerman, Drainage of curd, Ph.D. thesis, Landbouwuniversiteit Wageningen, Wageningen, The Netherlands, 1992. [10] J.C. Akkerman, F.H.J. Fox, P. Walstra, Drainage of curd: expression of single curd grains, Neth. Milk Dairy J. 48 (1994) 1–17. [11] J.C. Akkerman, F.H.J. Fox, P. Walstra, Drainage of curd: role of drainage equipment in relation to curd properties, Neth. Milk Dairy J. 50 (1996) 371–406. [12] O. Coussy, Mechanics of Porous Continua, Wiley, New York, 1995, p. 454. [13] R.W. Lewis, B.A. Schrefler, The Finite Element Method in the Deformation and Consolidation of Porous Media. Wiley, New York, 1998. [14] R. DeBoer, Theory of Porous Media: Highlights in Historical Development and Current State, Springer-Verlag, Berlin, 1999. [15] P. Dejmek, Personal communication, 1999. [16] D.S. Burnett, Finite Element Analysis, Addison-Wesley, Reading, MA, 1987. [17] C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, PA, 1995. [18] E. Tijskens, D. Roose, H. Ramon, J. De Baerdemaeker, Automatic differentiation for solving nonlinear partial differential equations: an efficient operator overloading approach, Numer. Algorith. 30 (2002) 259–301. [19] E. Tijskens, H. Ramon, J. De Baerdemaeker, FastDer++: efficient automatic differentiation for nonlinear PDE solvers, Math. Comput. Simul. 65 (2004) 177–190. [20] F. Mariette, Personal communication, 1998. [21] C. Lodaite, K. Östergren, M. Paulsson, P. Dejmek, One-dimensional syneresis of rennet-induced gels, Int. Dairy J. 10 (12) (2000) 829–834.