Bread dough rheology and recoil

Bread dough rheology and recoil

J. Non-Newtonian Fluid Mech. 143 (2007) 107–119 Bread dough rheology and recoil 2. Recoil and relaxation Roger I. Tanner ∗ , Shao-Cong Dai, Fuzhong Q...

2MB Sizes 1 Downloads 110 Views

J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

Bread dough rheology and recoil 2. Recoil and relaxation Roger I. Tanner ∗ , Shao-Cong Dai, Fuzhong Qi School of Aerospace, Mechanical and Mechatronic Engineering, University of Sydney, Sydney, NSW 2006, Australia Received 14 November 2006; received in revised form 1 February 2007; accepted 2 February 2007

Abstract It is shown that a simple Lodge model with a power-law memory function and a damage function which is a function of strain can describe the behaviour of a bread dough in steady shear, steady elongation, small and medium size sinusoidal strains, and shear stress relaxation. The number of parameters needed to be found from experiment is minimal, which is a clear advantage. Some questions remain about the frequency dependence in the sinusoidal oscillatory test. We also use the model to describe the difficult, but practically important problem of recoil after steady elongation. By a simple modification of the damage function concept for the recoil phase, we are able to describe the results of experiments at elongation rates between 0.001 and 0.1 s−1 and Hencky strains up to 2.5. Finally, although no explicit yield stress has been introduced into the model, results resemble those from models with a yield stress that depends on the rate of elongation. © 2007 Elsevier B.V. All rights reserved. Keywords: Bread dough; Viscoelasticity; Recoil; Relaxation; Damage function

1. Modelling of bread dough The present paper continues the exploration begun in Part 1 [1] of the use of a modified Lodge elastic fluid model with a damage function to describe the mechanics and rheology of bread dough. In Part 1 [1] we showed that the rheology of this soft, starch-filled solid could be described in suddenly started steady shear and elongation by a constitutive equation of the Lodge [2] form:  t σ + PI = f m(t − t  )C−1 (t  ) dt  (1) −∞

where ␴(t) is the stress at time t, P the pressure, I the unit tensor and C−1 (t ) is the Finger strain tensor [2,5] at time t computed relative to the configuration at the present time t. The memory function is assumed to be of the power-law form: m(t) = pG(1)t −(1+p)

(2)

where the constants p and G(1) can be found from small-strain oscillatory testing [1]. (Note that the constant G(1) here signifies the numerical value of the stress relaxation modulus G(t) at t = 1 s—the constant does not have the dimensions of stress.) The ∗

Corresponding author. Tel.: +61 2 9351 7153; fax: +61 2 9351 7060. E-mail address: [email protected] (R.I. Tanner).

0377-0257/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2007.02.001

“damage function” f is assumed to be a function of the Hencky strain at time t, computed relative to the state of rest, which is assumed to occur for times less than zero. For the JANZ dough used for the tests, we have already [1] reported on preparation methods, sinusoidal small-strain, shear and elongation beginning from rest, and on sinusoidal strain tests of moderate magnitude (∼10% strain). Here we present, for the same dough, relaxation after a shear strain applied at t = 0, and recoil in elongation from various initial strains. These experiments form, we believe, a unique set for a single dough, and we attempt to use Eqs. (1) and (2) to describe all the data. Two characteristics of the model are the economy of parameters, and the rapid decrease of the damage function f from unity at very small strains (∼0.001 or smaller) to a value of order 0.1 at strains of only 10–20% [1]. This remarkable softening with work input (kneading) is of course a well-known and valuable aspect of dough rheology. We begin by studying relaxation from a known shear strain, then consider the larger-strain oscillatory behaviour reported in Part 1 [1], and finally discuss recoil after elongational stress and release. 2. Stress relaxation Ideally, in this test a suddenly applied (step) of shear strain (γ 0 ) is applied at t = 0, and the decay of shear stress is measured.

108

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

Fig. 1. The shear strain–time curves achieved for various sized “steps” of shear. The maximum shear is used to normalize the curves. The ideal test would apply a true step at t = 0. The sizes of the steady-state shear steps (0.0001–0.3) are noted in the figure.

Fig. 2. Stress relaxation modulus after initial shear deformation. The final strain (γ 0 ) is only reached after t  1 s The decay curves are fairly well described by power-law curves of the form f(γ 0 )G(1)t−p for t  1 s. The ordinate is τ/γ 0 , the apparent shear relaxation function. The fitted lines all have the same slope t−0.27 .

In shear, Eqs. (1) and (2) can be written in the form:  t −p ˙  ) dt  τ=f G(1)(t − t  ) γ(t

relaxation tests (circles in Fig. 3). From Fig. 7 of Part 1 [1] we find the points shown as inverted triangles (steady shear) and squares (steady elongation). It is noticeable that for very long times (≥100 s) there is some softening present that is not predicted by the model. Looking at the small-strain oscillation results of Part 1 (Fig. 2) we see that these tests also show some deviation from a strict power-law at very low frequencies (<0.1 rad/s), and since measurements were not made below 0.06 rad/s, there may be effects at these low frequencies and long times which are not modelled correctly. However, for many process applications these long times/low frequencies will often be unimportant.We can also interpret these results using a KBKZ model. From well-known results [5] if the strain–time separation of the KBKZ kernel is assumed, then the

−∞

(3)

and the result for a step of shear of size (γ 0 ) is τ = γ0 f (γ0 )G(1)t −p

(4)

In experiments a true step of shear is not available, and Fig. 1 shows that around 1–2 s is needed to reach the final strain γ 0 in the Paar Physica MCR 300 instrument. There is also a variation of strain (and stress) with radius in the parallel plate apparatus, so here γ 0 and τ are actually taken as the strain and stress at 75% of the plate radius. It is known that this approximation is usually very little (O(1%)) in error [3]. If we idealize the response of Fig. 1 as a constant ramp γ˙ = γ0 /t0 , where t0 is the time to reach steady state and γ = γ 0 for t > t0 , we find the response at t (>t0 ) to be τ=

γ0 f (γ0 )G(1) 1−p [t − (t − t0 )1−p ] t0 (1 − p)

and for t  t0 we find:   t 2  p  t0  0 −p τ = γ0 f (γ0 )G(1)t 1+ +O 2 t t

(5)

(6)

Since in our case p = 0.27, the result (6) is within 1% of (4) for t/t0 ≥ 10. In fact the power-law asymptote will always be reached for large enough times. The present set of experiments shows behaviour close to the power-law relaxation at quite early times (Fig. 2), and in this figure we have fitted the curves shown of slope −0.27 to the stress data at t = 10 s. This enables us to find f(γ 0 ), since we assume G(1) = 12.2 kPa s0.27 , from the step strain of 0.01% (top curve in Fig. 2). (Previously we found G(1) ∼ 10.7 kPa s0.27 from oscillation; so there is reasonable concordance in these values.) Recalling that the Hencky strain εH ∼ (1/2)γ for small strains [4], we can plot f(εH ) from our

Fig. 3. Damage function f as a function of the Hencky strain εH for stress relaxation measurements ( ), finite-amplitude strain oscillations at 1 Hz (♦) and the average of 0.01, 1 and 30 Hz measurements ( ), steady shearing () and steady elongation results () from Part 1 [1] are also shown. For small εH (<0.05) the full line shows f = −(0.135 + 0.123 ln εH ); for larger εH values, the values for shear and elongation are shown in Fig. 7, Part 1 [1]; the formula for f in this range is given in Appendix A (Eq. (A4)).

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

result for a step strain size γ at t = 0 is τ = h(γ)G(t) γ

109

(7)

where h is the familiar damping function. Evidently both the KBKZ and damage function models give the same results for this deformation, and h = f, where f is given in Fig. 3. On the other hand, the predictions of models of the pom-pom variety [6,7] or the PTT model [5] in relaxation after step strain are known to be unrealistic [8] as the separable form shown in Eq. (7) does not occur; experiments [5] usually show the separable form of Eq. (7). Therefore, we shall not attempt to fit the data using these models, despite finding reasonable results for the recoil of polyethylene using the PTT models [9]. 3. Finite amplitude sinusoidal strain The damage function model predicts, for a shearing strain γ = γ 0 cos ωt, that G (γ 0 , ω) and G (γ 0 , ω) are related to the small-strain parameters G (ω), G (ω) by G (γ0 , ω) = fG (ω);

G (γ0 , ω) = fG (ω)

Fig. 4. Measured diameter at central plane of a specimen being elongated at a nominal rate of extension of 0.01 s−1 . The actual measured final rate of extension from the diameter data is about 0.0097 s−1 . For smaller strains (1 > εH ) the measured rate is about 0.0118 s−1 . The solid line shows the nominal strain rate of 0.01 s−1 .

(8)



and so if |G* | = (G 2 + G2 )1/2 , a similar rule applies to |G* |. Reference to Figs. 9 and 10 of Part I shows that Eq. (8) does not show the frequency dependence found in experiments. In Fig. 3 we have plotted the values of |G*(γ 0 , ω)|/|G* (ω)| taken at 1 Hz (ω = 6.28 rad/s) (♦ symbols) and the average over the values of the ratio at 0.01, 1 and 30 Hz ( symbols in Fig. 3). The value of Hencky strain to use in the damage function in this case is εH ∼ γ 0 /2, since εH is small here. The agreement between the values of f found from relaxation and the sinusoidal values is fair, although the latter are somewhat higher. We now consider the more complex problem of recoil. 4. Experimental data on elongational recoil In the recoil experiments, the specimens begin as discs (38 mm diameter × 8 mm thick) and a constant elongational strain rate is applied. To enable the stress to be determined, the specimen diameter at the mid-plane was measured using a video recorder. For the larger strains Fig. 4 shows that a reasonably close approximation to a uniform strain rate was achieved after a Hencky strain of order unity, and that the nominal rate of strain was within about three percent of the actual strain rate measured from the diameter. At smaller strains (εH ≤ 0.5) the error is larger, and the measured strain rate (Fig. 4) was about 19% larger than the nominal strain rate, due to the fairly complex transition of specimen shape at lower strains. At higher strains the specimens formed very satisfactory cylinders (Fig. 5). To alleviate the problem of distortion at low strains, specimens about 24 mm high and 14 mm in diameter were used for the εH = 0.2 and 0.5 recoil tests (Fig. 6). Because dough is comparatively insensitive to rate effects ((˙ε)0.27 , see Part I [1]) we have used the nominal strain rates (0.001, 0.01 and 0.1 s−1 ) in subsequent calculations. Errors of less than 5% in stress, even with smaller strains with 8 mm specimens are expected, and generally errors of around 1% are expected with larger strains.

Because of gravity and lack of symmetry in cutting, the recoil of top and bottom parts of the specimens were measured separately as a function of time and the total taken. (Variations of stress along the specimen due to gravity are only O(102 Pa), much less than the O(104 Pa) of the mean stress.) Fig. 7 illustrates the progress of a typical recoil. The recoverable strain εr is defined as   length after cutting (9) εr = ln length before cutting Three rates of elongation (0.001, 0.01, 0.1 s−1 ) and five strains at the cutting point (0.2, 0.5, 1.1, 1.81 and 2.5) were investigated and the recoil results are shown (symbols) in Figs. 8–12. There is a very rapid, nearly instantaneous recoil, followed by a very slow recovery. Because of the rapidity of recoil, we have considered wave propagation in the specimen. Following Pipkin [18], we find the propagation of a plane wave in a linear viscoelastic

Fig. 5. Specimen being elongated at a nominal strain rate of 0.01 s−1 at a Hencky strain of 2.

110

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

Fig. 6. For smaller strains a longer specimen was used. The act of cutting (1 s between pictures) is shown for the case ε˙ 0 = 0.1 s−1 at a Hencky strain of 0.5.

medium proceeds at a speed [ρJ]1/2 times a factor which is close to one. In the case of a longitudinal wave in a power-law material, we have the shear compliance J(t) = tp /G(1) so the longitudinal compliance is 1/3 of this. We calculate a wave speed U of [ρtp /3G(1)]−1/2 ·constant. In the case under discussion, ρ ∼ 1200 kg/m3 , p = 0.27 and G(1) ∼ 12 kPa s−p .

Hence U ∼ 5.5t−0.135 and for t = 0.001 s, U ∼ 14 m/s; for t = 1, U ∼ 5.5 m/s. Since the maximum strain before cutting is 2.5, the length of each specimen half is about 50 mm maximum, so any wave propagation effects are over in less than about 5 ms, which cannot be resolved in the experiments. Hence, we have neglected wave propagation in the specimens. We now turn to the analysis of the data.

Fig. 7. A typical sequence showing the cutting of the specimen at ε0 = 1.81, ε˙ 0 = 0.001 s−1 . The pictures show recoil at 0, 1, 5, 10, 30, 60 and 180 s after cutting.

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

Fig. 8. Recovery curves for ε˙ 0 = 0.001, 0.01 and 0.1 s−1 stretching to ε0 = ε˙ 0 t = 2.5. The recovery curve uses Eq. (50) for the calculation. The formula for f is given in Eq. (A4); here f ∼ 0.0327. Note the initially rapid recovery and then the slow recovery. The shapes of the recoil curves are well described by the calculations.

111

Fig. 10. As in Fig. 8, except initial strains are 1.1. f = 0.0672 here.

Fig. 11. As in Fig. 8, except initial strains are 0.5. f = 0.0874 here.

Fig. 9. As in Fig. 8, except for initial strains which are 1.81, 1.81 and 1.71 for ε˙ 0 = 0.001, 0.01 and 0.1 s−1 , respectively.

5. Recoil calculations—linear cases In Part 1 [1] we showed three models that have been proposed for dough: (a) Lerchental and Muller’s (1967) model is extremely complex and it is hard to see how to write down a proper threedimensional generalization; it contains springs, dashpots, yield elements and shear-pin elements. (b) Bloksma’s (1971) model is much simpler and contains only one yield element. However, even for small strains it contains only two relaxation times, which are clearly insufficient to describe recoil. (c) The (2000) model of Phan-Thien et al. [10] consists of a nonlinear spring in parallel with a number of PTT elements.

Fig. 12. As in Fig. 8, except initial strains are 0.2, 0.2 and 0.23. f = 0.1139 for ε0 = 0.2 and f = 0.1089 for ε0 = 0.23 here. The lowest curve shows almost perfect recoil.

112

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

Recovery/recoil is always 100% with this model. Hence, there is no convincing treatment of dough recoil available despite its importance. In the following, we explore the use of the damage function and KBKZ models in recoil and thereby we avoid the explicit use of yield elements. Recoil presents a difficult situation to analyse and before considering the analysis of the experiments above, we give some illustrative background examples. This section is intended to demonstrate the effects of the damage function and the power-law memory function on recoil in a simple way. The points demonstrated are: (a) the damage function can limit recoil even with a linear spring; (b) a power-law kernel always gives complete recovery in the linear viscoelastic case; (c) an initial jump of strain occurs with a Maxwell model; (d) there is a similarity in behaviour between Wagner’s irreversible model and the damage function model. 5.1. Elastic model To begin, let us consider a simple elastic spring (modulus E) subject to damage. The constitutive model is assumed to be: ␴ = fE␧

(10)

Let us suppose the material is composed of many filaments and a proportion of these break at strain ε. Hence, we assume f = f(ε) in Eq. (10). At the beginning of recoil, the broken strands are unstressed. Let the maximum strain be εmax (see Fig. 13). At the end of extension: ␴ = f (␧max )E␧max

(11)

Upon recoil, the fraction (1 − f) of damaged strands is compressed; the total compressive stress (negative) in these strands is ␴c = (1 − f )Er ␧r

(12)

where ␧max > ␧r . The recoil modulus Er could be different from E, the loading modulus. The total stress after recoil is zero, so σ = fE(εmax − εr ) − (1 − f )Er εr = 0

(13)

Hence we find: εr Ef = εmax Ef + Er (1 − f )

(14)

Note that when f = 1 (no damage) εr = εmax (complete recoil); when f = 0 (fracture) εr = 0 (no recoil). If Er = E, then εr = fεmax . 5.2. Linear viscoelastic Maxwell model Let



␴ = 3f (␧)

t

−∞

˙  )(t  ) dt  G(t − t  )␧(t

(15)

Suppose steady extension occurs at a rate ε˙ 0 for t = 0 to t0 , and recoil occurs after t0 . The stress at t0 is then  t0 G(t − t  ) dt  (16) σ(t0 ) = 3f (˙εt0 )˙ε0 0

If we imagine that the damaged fraction (1 − f) is unloaded at this point, but upon retraction begins to bear a compressive stress, then the total stress for t > t0 is  t  t σ(t)=3f G(t − t  )˙ε(t  ) dt  + 3(1 − f ) G(t − t  )˙ε(t  ) dt  t0

0

where f = f (˙ε0 t0 ). Hence, for t > t0 :  t0  t 0 = f ε˙ 0 G(t − t  ) dt  + G(t − t  )˙ε(t  ) dt  0

t0

If we assume a single-mode Maxwell model, where  t η G(t) = exp − λ λ

(17)

(18)

(19)

and an instantaneous unloading at t = 0 [9], with a strain jump εr , then from (18) we find: εr = f (˙ε0 t0 )(λ˙ε0 ) {1 − e−t0 /λ }

Fig. 13. Purely elastic spring with a damage function showing recoil.

(20)

A sketch of the situation is shown in Fig. 14. For multiple (discrete) relaxation times there will be an initial jump and then a slow recovery [9], which is more complex to analyse. In the linear Maxwell case, for a rapid loading, where t0 /λ 1, εr = f ε˙ 0 t0 (elastic limit), and for very slow loading εr ∼ 0 (viscous limit). Here there is always an imperfect recoil. In Eq. (17) we have assumed that the constitutive law for recoil is the same as that for loading, that the entire (1 − f) fraction of the material participates in the compression loading in the recoil process. In fact, in a later section we will assume that only a multiple β of the (1 − f) damaged strands participates in recoil (0 ≤ β ≤ 1).

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

zero:



τ(t, γ) =

t

−∞

G(t − t  )

113

d (γf ) dt  dt 

(27)

which is actually linear in the strain measure χ = γf. Clearly, after any shear stress history imposed in the range 0–t0 we can follow our previous reasoning to show that for a power-law memory function, γf must always return to its original state (γ = 0) and hence recoil is always complete. Whilst this problem can be avoided by using Wagner’s irreversible kernel [5,11] this is really very close to using a damage function approach; actually we could combine the two to give  t τ(t, γ) = f m(t − t  )γh dt  (28) −∞

where f is “irreversible”, and h(γ) is reversible, but we shall not do that here. If we use the KBKZ approach, it also follows that elongational recoil will always be perfect. We now consider the use of Wagner’s irreversible kernel idea [11].

Fig. 14. Sketch of recoil of a linear viscoelastic Maxwell model.

5.3. Linear power-law kernel—no damage function In this case, G(t) = G(1)t−p , and we have  t σ = G(t − t  )˙ε(t  ) dt  3 0

5.4. Wagner irreversible kernel (21)

By introducing the creep function J(t), we can compute the strain ε as [5]:  1 t dσ ε(t) = J(t − t  )  dt  (22) 3 0 dt Now if σ = 0 for t > t0 , then  1 t0 dσ ε(t) = J(t − t  )  dt  3 0 dt Integration by parts gives, since σ(t0 ) = σ(0) = 0:  1 t0  σJ (t − t  ) dt  ε(t) = 3 0 For large t  t0 , J (t − t ) ∼ J (t), and so for large t:  t0 1 σ dt  ε(t) ∼ J  (t) 3 t

(23)

(24)

For moderate strains (ε < 0.1) we will use, for simplicity, a small-strain form in shear:  t ˙  ) dt  τ= G(t − t  )hγ(t (29) −∞

where the damping function h is a function of the strain γ. In the “irreversible” calculation [5,11], during recoil h remains at its minimum value, which occurs when t = t0 and t = 0. Eq. (29) is thus nonlinear only through the h-function [12]. For a steady shearing of γ˙ 0 imposed at t = 0 and ceasing at t = t0 , the recoil equation for shear stress is  t0  t ˙  ) dt  G(t − t  )h dt  + hm G(t − t  )γ(t (30) 0 = γ˙ 0 0

t0

where hm is the minimum value of h, to be used during recoil, and in the first integral h = h(γ˙ 0 (t − t  )). If G(t) has the form (2), and h = f, from Fig. 3 (noting εH = γ/2 here): (25)

In the case G(t) = G(1)t−p , J (t) = tp−1 /(p − 1)!(−p)!; and since p < 1, ε(∞) → 0, or recovery is complete. Thus, the power-law material with no damage (f = 1) behaves elastically (εr = ε˙ t0 ); when f < 1 we expect the recoil to be governed by the elastic law (14), so that if elongation and recoil have the same moduli G(1), then εr is expected to be roughly f ε˙ 0 t0 . We can also investigate recoil for the KBKZ model. We have, for shear [5]:  t τ(t, γ) = m(t − t  )f (γ)γ(t  ) dt  (26) −∞

similar to Eq. (1); here γ is the strain at t relative to the present time t. Integration by parts gives, noting that γ(t) and G(∞) are

h = A + B ln γ

(31)

then the first integral in (30) can be evaluated. We find A = −0.042, B = −0.113, and   t p 1−p  γ˙ 0 γ0 B −p ˙  ) dt  (32) hm − + hm (t − t  ) γ(t 0= 1−p 1−p t0 where hm is the minimum h, equal to A + B ln γ 0 . Eq. (32) is clearly similar to the damage function formulation (18). The correspondence is exact if we set f =1−

B hm (1 − p)

(33)

Thus if γ 0 = 0.1, the equivalent f = 1.71. However, in the damage model f cannot exceed 1, so the result (33) is unwelcome, and leads to a recoil greater than that for the damage function model.

114

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

This is disappointing, as the irreversible model works well with polyethylene [5,11]. The greater elasticity of dough at low rates of deformation is probably the cause of this. Hence, we will return to the damage function model. We note, from these examples, that the damage function approach can yield an imperfect recoil even in highly elastic systems. 6. Nonlinear recoil calculations—an approximate treatment Generally, the treatment of the full nonlinear response will need numerical analysis, but before doing that we present a simple treatment of Eqs. (1) and (2). There are considerable difficulties in solving the exact nonlinear recoil problems, but we note that (Fig. 3, Part 1 [1]) the material is almost elastic for large elongational strains, and hence the type of simple model discussed in Section 5.1 may be useful. Let us suppose the specimen is stretched at a uniform rate ε˙ 0 for the period t = 0 to t0 , so that the Hencky strain before recoil ε0 = ε˙ 0 t0 . If the specimen length was at rest, then at t = t0 it will have a length exp(ε0 ). Recoil then takes place (Fig. 8), initially very rapidly, and then slowly, over a long period (O(103 s)). The damage function f(ε0 ) enables the stress σ 0 to be found [1]: σ0 = pG(1)(˙ε0 )p F (ε0 )f (ε0 )

(34)

where the function F is given in Fig. 3 of Part 1 [1], and f is given in Fig. 3. Upon recoil, as in Eq. (14) above, we suppose f = 1. If the material were a rubber with damage, then reloading from an unknown strain εa so as to reach σ 0 when ε = ε0 in the undamaged material (f = 1), would give a method of finding εa , the residual strain. The measurements give the recoverable strain εr as in Eq. (9). Now if the specimen is originally of length , then after extension, the length is exp(ε0 ), and (Fig. 15) after recoil it is of length exp(εa ). Hence, the recoil strain εr is εa – ε0 , a negative quantity. On reloading, the specimen has a length exp(εa ) and is stretched to exp(ε0 ), so the specimen strain in this case is ε0 – εa , the negative of the recoil strain. On the unloading path we have no control over the rate of extension, which actually varies continuously. For a simple model, we will suppose that recoil takes place at a constant average strain rate ␧¯˙ which is not known. Then using (34) we find: p

F (εr )(ε¯˙ ) = (˙ε0 )p F (ε0 )f (ε0 )

(35)

Fig. 15. “Elastic” model of recoil. Recoil strain εr = εa − ε0 ; ε0 = ε˙ 0 t0 . The process of recoil is replaced by a reloading at a constant average rate from εa to ε0 .

Looking at the experimental recoil data at early times we see that (Figs. 8–12): (i) when ε˙ 0 is larger, so is the recovery rate; (ii) when ε0 is larger, so is the recovery rate and at later times recovery is very slow. Neither of these conclusions is unexpected. A reasonable empirical description of the data is given if we set (˙ε in s−1 ):  p 7.0˙ε0.16 ε˙ 0 0 = (36) ε0 ε¯˙ Calculations showing the agreement with experiment for selected extensions appear in Table 1. For a complete lack of recoil, εr = 0; for a complete recoil, εr = −ε0 . While an empirical fit to the total recoil can be described by this method in the range of tests given in Table 1, the strain–time behaviour is not correctly modelled, and hence we now examine the recoil more closely. 7. Nonlinear calculations of jump strain It is noticeable in the experimental data presented in Figs. 8–12 that there is an appreciable sudden jump strain or recoil, especially at higher rates of extension. We may analyse

Table 1 Comparison of measured and computed recoil ε˙ 0 (s−1 )

ε0 2.5

0.001 0.01 0.1

1.8

0.5

Experiment

Eq. (35)

Experiment

Eq. (35)

Experiment

Eq. (35)

0.45 0.68 0.86

0.52 0.70 0.90

0.44 0.50 0.66

0.3 0.43 0.62

0.1 0.17 0.25

0.1 0.24 0.36

Figures given are for |εr | (εr is actually negative).

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

115

Table 2 Jump strains (λm = 0.001 s, p = 0.27) ˙ ≡ ␧H Hencky strain ␧t

Rate of extension (s−1 )

2.5

0.0001 0.01 0.1

1.81

0.5

εj (experiment)

εj (calculated)

εr

εj (experiment)

εj (calculated)

εr

εj (experiment)

εj (calculated)

εr

−0.26 −0.49 −0.67

−0.46 −0.63 −0.83

−0.45 −0.68 −0.86

−0.19 −0.33 −0.47

−0.20 −0.32 −0.47

−0.44 −0.50 −0.59

−0.02 −0.06 −0.17

−0.026 −0.047 −0.085

−0.1 −0.17 −0.26

Note: All strains are actually negative. The 0.1 s−1 results for the medium strain refer to εH = 1.71, not 1.81 (calculations always used the same strain as the corresponding experiments).

˙ The and a similar formula for τ 1 , if 2␧˙ is substituted for −␧. value of λm will be small—smaller than the inverse of the maximum frequency reached in our oscillatory testing (628 rad/s), and hence we will assume ␭m ∼ O(10−3 s). To evaluate the integral in (41) it is convenient to let t − t = s, and expand the exponential in a power series in ε˙ s and integrate term by term, finding: 

 p˙εt p˙ε2 t 2 p˙ε3 t 3 −p −p −˙εt 4 τ2 = G(1) λm − t 1−e + (42) − + − O(t ) 1 − p 2(2 − p) 6(3 − p)

this jump by following Lodge [2], see also [9]. Consider Eq. (1) and assume the strain history is a uniform extension rate (˙ε0 ) for t > 0, and then a sudden recoil at time t immediately after the dough is unloaded. Generally, a much slower recovery will also follow (Fig. 8). Let the jump strain in the axial direction be of magnitude ε (so that the specimen strain along the axis after

the jump is ε times the maximum strain ε˙ 0 t). Then if we denote (from Eq. (1)) ␶1 = ␴11 + P, ␶2 = ␴22 + P, where the elongation is in the one-direction, we find that the jump strain (ε) will unload the specimen if [9]:  1/3 τ2 ε= (37) τ1 where the stresses are evaluated just before the jump. In the present model:  t −1  −1  [τ1 ; τ2 ] = f m(t − t  )[C11 (t ); C22 (t )] dt  (38) −∞

and 

−1 C11 = e2˙ε0 (t−t ) ,



−1 C22 = e−˙ε0 (t−t )

(39)

Hence, f does not come into the calculation of the jump strain. However, using the kernel (2), it is easy to see that the integrals for τ 1 and τ 2 do not converge, and the calculation cannot be made. To avoid this problem, we note that the relaxation spectrum corresponding to (2) is H(λ) =

G(1)λ−p (p − 1)!

(40)

In this section we will use a modified relaxation spectrum and assume that there are no modes shorter than a minimum λm in (40); equivalently, we set G(t) = G(1)t−p for t > λm , and G(t) = −p G(1)λm for 0 ≤ t ≤ λm . Because the m-function is the derivative of G, m = 0 for 0 ≤ t ≤ λm , and equal to Eq. (2) for larger times. We now find finite values for τ 1 and τ 2 , and the calculation can proceed. We find:  t−λm −(p+1) ε˙ t   −˙εt −p −˙εt (t − t  ) e dt τ2 = G(1) e t + pG(1) e 0

(41)

and a similar formula for τ 1 found by substituting 2␧˙ for −␧˙ (for the experimental conditions investigated τ 1 and τ 2 are given by the above formulas to about 1% accuracy); λm is not known; we assume a value of 0.001 s; p = 0.27 as before; the exact value of G(1) is not important. Note that Eq. (42) is a convergent series for all values of ε˙ t. We can now compute the jump strain from (37); to compare with experiments we use the Hencky strain:   1 τ2 ln ε = εj = ln (43) 3 τ1 To compare with experiment we need to measure the jump strain from Figs. 8–12. To do this the jump strain is measured from the cutting moment t = 0 to the first dot that departs from the vertical axis in Figs. 8–12—a time of around 5 s. Let εH = ε˙ 0 t. Then, following Table 1 we have (Table 2) to tabulate both jump (␧˙ j ) and total (εr ) recoils for the measurements made. Clearly the slowly loaded specimens have a lesser jump and a larger slow recoil. Given the various theoretical and experimental difficulties encountered, the agreement is encouraging. For very large strains, at finite strain rates, the following asymptotic formula can be found:   1 2 t ln ε ∼ p ln − εH (44) 3 λm 3 so that ultimately the recoil is 2/3 of the initial strain, and there is always a slow recoil, back to εH . If we permit λm to become smaller and smaller, from (42) we see that ultimately τ 1 → τ 2 and the jump strain disappears. 8. Numerical calculations of recoil To compute the complete recoil, we need to use a numerical method. Although it is possible to compute using the integral

116

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119 Table 3 Discrete 16-mode and 14-mode spectra for dough

Fig. 16. G and G comparing measured and 14-mode and 16-mode results: (circles) experiments; (solid lines) 14 modes and (dashed lines) 16 modes.

λi (s)

gi (Pa)

0.000316 0.001 0.00316 0.01 0.0316 0.1 0.316 1.0 3.16 10.0 31.6 100.0 316.2 1000.0 3162.0 10000.0

32421.0 23759.0 17411.93 12757.41 9350.76 6851.13 5021.65 3679.28 2696.79 1975.89 1448.26 1061.11 777.63 569.85 417.61 1140.0

G(1) = 10.7 kPa, p = 0.27.

form (1), it is simpler to use a differential form of (1), using a set of discrete-mode Maxwell models. However, to be sure that the discrete-mode model mimics the integral form (1), we have computed the response in small-strain oscillation and in elongation in Figs. 16 and 17. Fourteen modes had to be used to obtain a satisfactory response; 16 modes were better. Relaxation times were chosen at half-decade intervals from 104 to 0.003162 s for 14-mode (or 0.0003162 s for 16-mode), and Table 3 shows the data used. Baumg¨artel and Winter [13] have shown how to make a discrete spectrum from a continuous one. If the discrete spectrum result for G(t) is G(t) =



gn e−t/λn

(45)

n

then with a constant spacing of half a decade we have [13]: gn =

1 H(λn ) ln 2



λn−1 λn+1

 (46)

Fig. 17. Linear relaxation functions of 14-mode and 16-mode Maxwell models, showing they agree with the power-law for t > 0.001 s.

where the H-values are given by Eq. (40). A slightly better rule with the power-law spectrum (40) finds: gn = H(λn )

r p/2 − r −p/2 p

(47)

where r = λn+1 /λn . The advantage of (47) is that the longest relaxation time can have all the results out to infinity lumped with it; in the present case r = 3.162, p = 0.27, the final gn (λ = 104 s) has 3.742 times the weighting expected from rules (46) and (47) (Table 3). The results for the 14-mode and 16-mode G , G response are shown in Fig. 16, which was compared with oscillatory data out to 100 Hz, and the corresponding G(t) are shown in Fig. 17. Note the power-law behaviour is not seen below about 10−3 s. The calculated elongation stresses using a discrete spectrum (no damage, f = 1) are also compared to those found using a continuous spectrum (Fig. 18). The computations of recoil use the damage function concept as discussed above, and the methods of Ref. [9] described in Appendix A; 14-mode

Fig. 18. Comparing the integral model (Eqs. (1) and (2)) and 14-mode upperconvected Maxwell model in elongation. Points: 14-mode UCM; solid lines: integral model (Eqs. (1) and (2)).

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

and 16-mode results for recoil were indistinguishable, and we report only 14-mode results. If the damage function f = 1, then the recoil is almost perfect. As we have seen above, we need to consider what the material does after cutting and during the recoil. The simplest assumption is that the (1 − f) part of the material behaves like the undamaged material, but that it is unstressed just before and after cutting; this assumption was used in Section 5.2 (Eq. (17)) for the linear case. The recoil with this assumption is predicted to be smaller than observation for ε˙ 0 t0 = 1.8 and 2.5 when the f-values of Fig. 3 are used. For ε˙ 0 t0 = 0.5 the measured recoil is again larger than the predicted recoil, and this is also the case for smaller strains, where often an almost complete recoil has been reported [14]. To address these problems we return to Bueche’s original Mullins effect concept [15]. He assumed that as elongation of various strands proceeded, that breakage occurred, and upon retraction the broken chains do not participate in generating stress. If we consider Eq. (1) in the case of 1 − C −1 ), by S(t ) we have, during elongation, then denoting (C11 22 loading:  t σ = f (ε) m(t − t  )S(t  ) dt  (48)

117

Fig. 19. The function β for various ε˙ 0 and ε0 . The lines are given by Eq. (50).

−∞

For the recoil, where σ = 0 and t > t0 , the stress is composed of two parts: (a) The component due to unbroken strands, equal to t f (εmax ) −∞ m(t − t  )S(t  ) dt  where εmax is the maximum strain. (b) The “broken” and new strand contribution on recoil. t The second component will be (1 − f ) −∞ m(t − t  )S(t  ) dt  if we assume that all (1 − f) strands contribute fully to recoil, as assumed above. Similar to the elastic argument given in Section 5.1, we will use a factor β to multiply (1 − f) to describe this contribution to the stress. Equating the total stress to zero, we find:  t0  t 0=f mS dt  + [β + (1 − β)f ] mS dt  −∞

+ β(1 − f )s(t0 )



t0

t0

−∞

m dt 

(49)

If β = 0 (no participation in the recoil stress by the 1 − f portion), it is easy to see that full recovery takes place; if β = 1, then results with a finite recoil can be found. As mentioned above, for ε0 = 0.5 the predicted recoil is too small. We have suggested that β is a modifier of the strand (1 − f) fraction, but this is not the only possibility. For example, a survey by Muhr [16], discussing the Mullins effect in filled rubbers, suggests that the relaxation function in recoil is different from that in extension. Hence, one could also regard β as a modifier of G(t) or m(t) in the recoil phase or, we could also assume that it reflects a balance between reversible and irreversible events in the network. In the present work we shall regard β as a (positive) function of strain and strain rate. It is known [14] that for εmax ≤ 0.25 and quick loading the recoil is observed to be almost perfect, so we see that β must be near zero up to this strain. Attempts to use β-

Fig. 20. Recovery strain predictions as a function ε˙ 0 and ε0 . Points: experiments; lines: predictions.

values which depend only on the Hencky strain ε0 (at t = t0 ) show insufficient change with strain rate ε˙ 0 ; hence we are forced to assume that β is a function of both ε0 and ε˙ 0 . A simple form for β (Fig. 19) is β=

1 35 (2ε0

+ 1.65)(1 − log ε˙ 0 ),

ε0 ≥ 0.5,

β = 0.45ε0 − 0.0757 log ε˙ 0 − 0.1493,

ε0 < 0.5

(50)

Note: Whenever Eq. (50) is negative, β must be set to zero. With Eq. (50), the results are well-described (Figs. 20 and 8–12 full lines), including the strain–time recovery histories. 9. Conclusions We have proposed a simple model (Eqs. (1) and (2)) of bread dough behaviour and have tested it in many deformation patterns. The idea of a damage function f, depending on the maximum strain reached relative to the rest state is introduced. The f-concept permits us to describe steady shear and elongation, as described in Part 1 [1] by using relatively few parameters. The

118

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

parameters needed are: (a) The slope parameter p and the value of G(1), both of which can be found by sinusoidal testing at very small strains (O(0.001)). (b) The f-function. A single function may be found to describe shear and elongation [1] here we have presented the data (Fig. 3) as a function of the Hencky strain εH , with the same curve for elongation and shear. A convenient method of finding f for small strain is via shear stress relaxation tests (Section 2). For large strains up to fracture, elongation and/or shear tests are required. These tests can be made in many rheological laboratories – however some additional data for recoil description are also needed – see below. Work on various dough mixes is needed to establish how the rheological parameters vary with wheat type and dough composition. Using this model, successes include: (1) Steady shear and elongation descriptions up to the fracture point. (2) The prediction of the relaxation of shear at modest shear strains (up to 30%). (3) The correct “instantaneous” jump behaviour after stress release in elongation. Here it is necessary to limit the shortest relaxation times used to 0.001 s, to avoid divergent integrals. The prediction of recoil is more difficult (Fig. 20). If we assume that the (1 − f) fraction of the material is unstressed after stress release, due to damage, and upon recoil begins to undergo stress beginning from zero stress at the time of the start of recoil, then the predicted total recoil is generally too small, and the effect of the rate of loading is not correctly reflected in the recoil. Unfortunately, we have not found any useful guide to these effects in the literature [16]. In order to address these concerns, we assume that a factor β(˙ε0 , ε0 ) (Fig. 19) multiplies the (1 − f) factor for the recoil; the rest of the (undamaged) material is not altered. We find that: (a) β is a decreasing function of ε˙ 0 ; (b) β is an increasing function of ε0 . With regard to (a), we suggest that larger initial strain rates leave the damaged strands with less time to re-establish bonds with starch particles, and so there is a general softening of the memory function. Although we have represented this by a factor 1 − log ε˙ 0 in Eq. (50), an almost equally good description is given by a multiple of ε˙ −0.17 . In any case, for very large strain 0 rates, β tends to zero, and recoil becomes more complete. With regard to (b), we suggest that the larger strains occurring before damage enable the formation of entanglements, thus, increasing β and decreasing recoil. Obviously, a proper micromechanical theory of these effects is needed. If the damage function f is set equal to 1 (no damage) then the recoil is total or very nearly so. It was already noted (1932–1937) by Schofield and Scott Blair [14] in pioneering

papers on dough rheology that for elongational strains of less than 30% (εH ∼ 0.24) recovery appeared to be complete. For the theory proposed here, we see that for small strains (Section 5.3) recovery is expected to be complete; if we use the ideas of Section 6, for slow, small strains, then complete recovery again seems possible. These conclusions are in rough agreement with Schofield and Scott Blair’s observations. We can now look back at the finite-amplitude oscillation results (Section 3). It seems that this kind of test may give us an insight into the modification of the memory function in reversing strains, and hence lead to a resolution of the frequencydependence of the softening, which is not captured by the damage model. Again, no useful theory seems to be at hand. Although we used a discrete spectrum for the recoil calculations, the discrete spectrum can be deduced from Eq. (2) by lumping the continuous spectrum at the spectral lines, so that no extra information or numerical analysis is needed to find the Maxwell spectrum. It has not been necessary to introduce a yield stress (see [17]), although in some regards the damage function acts in a similar way as an irreversible element. Looking at Fig. 20, we see that recoil is nearly complete for a loading rate of 0.1 s−1 and a total strain of 0.23. This corresponds to a maximum stress of 1 kPa, which is higher than the yield stress quoted by Bloksma [19]. On the other hand, for a slow loading (˙ε0 = 0.001 s−1 ) the maximum stress is only about 300 Pa, but some permanent deformation ensues; thus, the “yield stress” is lower than the figure given in [19]. While it is clearly possible to have a rate-dependent yield stress, it introduces complexity. Here we have not needed the yield concept. Yield-like behaviour arises automatically. Acknowledgements We thank the Australian Research Council for a Discovery Grant which enabled the work to be done. We also thank Dr. Ferenc Bekes (CSIRO and FBD Ltd.) for his assistance in selecting the dough and advice on dough mixers and mixing. Appendix A. Calculations of recovery strain In present paper, the recovery strain after elongation for bread dough has been calculated. In our investigation, the UCM model has been used as the constitutive equation with the help of a damage function f and the function β multiplying the fraction of damaged strands (1 − f). For any single mode (λi , gi ), the equation of the UCM model is given by λi

␶i + ␶i = 2ηi d Δt

(A1)

where / t is the usual upper-convected time derivative [5], ␶i is the stress tensor due to the ith

mode and the total stress is the sum of the stress ␶i , (␶ = i ␶i ), λi is a set of relaxation times, ηi is the corresponding viscosity, and d is the rate of deformation tensor equal to (L + LT )/2. The velocity gradient tensor L is defined as L = (v)T or Lij = ∂vi /∂xj , where v(x, t) is the velocity field at place x at time t. The upper-convected

R.I. Tanner et al. / J. Non-Newtonian Fluid Mech. 143 (2007) 107–119

derivative / t is defined as

␶ ∂␶ + (v∇)␶ − ␶LT − L␶ ≡ ∂t

t

and the total strain (relative to the rest state at t = 0) is given by (A2)

We used 14 modes in the calculations; increasing this to 16 modes gave results practically identical to the 14-mode case. The complete calculations of recovery strain of bread dough after elongation have been carried on two steps. • Step 1. Stretch (loading): The specimen of dough is stretched with a strain rate ε˙ 0 from time t = 0 to t0 starting from rest. For this process, the stress state can be obtained by numerically solving Eq. (A1) with the damage function f at time t (0 ≤ t ≤ t0 ). The total stress is given by UCM UCM σ = σ11 − σ22 = f (σ11 − σ22 ) = fσ UCM

(A3)

where σ UCM denote the result obtained by solving Eq. (A1) and f is the damage function. We have assumed that f is only a function of the Hencky strain ε = ε˙ t and it has the form [1]:

1 30 f= − 0.0225ε + 0.087 e−(0.32ε) 1.1026 1.0953 + (111.739ε) (A4) • Step 2. Recoil (after unloading): At the time t = t0 , we suddenly remove the load and then recoil is started. The recovery strain (Hencky strain) can be calculated by numerically solving Eq. (A1) with the damage function f and function β, starting from t = t0 . We have found that β not only depends on the total Hencky strain ε0 at the time t = t0 , but also depends on the stretch strain rate ε˙ 0 . A simple form of β is shown as Eq. (50). Note both f and β are constant for the recoil process (t > t0 ). For any time step interval tN , we assume that the deformation rate ε˙ N is a constant, here N is the number of the iteration. Therefore, the time step intervals t have to be very small, especially for the early stages of recoil, and t = 10−6 (s) has been used in our calculations. By numerically solving Eq. (A1), to find a critical deformation rate ε˙ N (negative) to satisfy the condition of the total stress: UCM σ = fσfUCM + β(1 − f )σ1−f =0

(A5)

where σfUCM is obtained by solving Eq. (A1) with the initial stress at t = t0 equal to the stress σ UCM just before unloading UCM is the result by solving Eq. (A1) but the iniat t = t0 , σ1−f tial stress at t = t0 is set to be zero since the partial strands are assumed to start compression after unloading at t = t0 . Once Eq. (A5) is satisfied, the critical deformation rate ε˙ N for the Nth time interval is found. The increment of the recovery strain is therefore obtained ( εr )N = tN ε˙ N . The total recovery strain is the sum of strain increments from t = t0 : εr =

N j=1

( εr )j

119

(A6)

ε = ε0 + ε r

(A7)

Note the recovery strain εr is negative. References [1] R.I. Tanner, F. Qi, S.-C. Dai, Bread dough rheology and recoil. 1. Rheology, J. Non-Newt. Fluid Mech, in press. [2] A.S. Lodge, Elastic Liquids, Academic Press, London, 1964. [3] M.T. Shaw, Z.Z. Liu, Single-point determination of nonlinear rheological data from parallel-plate torsional flow, Appl. Rheol. 16 (2006) 70–79. [4] V. Kitoko, M. Keentok, R.I. Tanner, Study of shear and elongational flow of solidifying polypropylene melt for low deformation rates, Korea–Aust. Rheol. J. 15 (2003) 63–73. [5] R.I. Tanner, Engineering Rheology, 2nd ed., Oxford University Press, 2000. [6] R.J. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag–strain coupling in branched polymer melts, J. Rheol. 44 (2000) 121–136. [7] W.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Viscoelastic analysis of complex polymer melt flows using the extended pom-pom model, J. NonNewt. Fluid Mech. 108 (2002) 301–326. [8] R.I. Tanner, On the congruence of some network and pom-pom models, Korea-Aust. Rheol. J. 18 (2006) 9–14. [9] R.I. Tanner, A.M. Zdilar, S. Nasseri, Recoil from elongation using general network models, Rheol. Acta 44 (2005) 513–520. [10] N. Phan-Thien, M. Newberry, R.I. Tanner, Non-linear oscillatory flow of a soft solid-like viscoelastic material, J. Non-Newt. Fluid Mech. 92 (2000) 67–80. [11] M.H. Wagner, Elongational behaviour of polymer melts in constant elongation-rate, constant tensile stress, and constant tensile force experiments, Rheol. Acta 18 (1979) 681–692. [12] A.D. Drozdov, Mechanics of Viscoelastic Solids, Wiley, Chichester, 1998. [13] N. Baumg¨artel, H.H. Winter, Interrelation between discrete and continuous relaxation spectra, in: P. Moldenaers, R. Keunings (Eds.), Theoretical and Applied Rheology, Proceedings of the XI International Congress Rheology, Brussels, 1992, pp. 893–895. [14] R.K. Schofield, G.W. Scott Blair, The relationship between viscosity, elasticity and plastic strength of soft materials as illustrated by some mechanical properties of flour doughs. I, Proc. Roy. Soc. London A 138 (1932) 707–719; R.K. Schofield, G.W. Scott Blair, The relationship between viscosity, elasticity and plastic strength of soft materials as illustrated by some mechanical properties of flour doughs. II, Proc. Roy. Soc. London A 139 (1933) 557–566; R.K. Schofield, G.W. Scott Blair, The relationship between viscosity, elasticity and plastic strength of soft materials as illustrated by some mechanical properties of flour doughs. III, Proc. Roy. Soc. London A 141 (1933) 72– 85; R.K. Schofield, G.W. Scott Blair, The relationship between viscosity, elasticity and plastic strength of soft materials as illustrated by some mechanical properties of flour doughs. IV, Proc. Roy. Soc. London A 160 (1937) 87– 94. [15] F. Bueche, Molecular basis for the Mullins effect, J. Appl. Polym. Sci. 4 (1960) 107–114. [16] A.H. Muhr, Modeling the stress–strain behavior of rubber, Rubber Chem. Technol. 78 (2005) 391–425. [17] G.E. Hibberd, N.S. Parker, Nonlinear creep and creep recovery of wheatflour doughs, Cereal Chem. 56 (1979) 232–236. [18] A.C. Pipkin, Lectures on Viscoelasticity Theory, 2nd ed., Springer-Verlag, New York, 1986. [19] A.H. Bloksma, Dough structure, dough rheology, and baking quality, Cereal Foods World 35 (1990) 144–237.