The fluid dynamics of the chocolate fountain

We attempt to solve the Navier–Stokes equations, then, equation (1), along with mass conservation for an incompressible fluid,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot {\bf{u}}=0.\end{eqnarray} \tag{ 8 }$$

The radius of the pipe is a, and the natural choice is to use cylindrical coordinates (r, θ, z), so ${\bf{u}}=({u}_{r},{u}_{\theta },{u}_{z}),$ etc.

Symmetry in the problem leads us to expect no motion in the r- or θ-directions, and to expect axial velocity to be merely a function of radial position, u_z = u_z(r). Furthermore, searching for steady flow solutions sets all time derivatives to zero. Under these conditions, the only non-zero components of shear rate are

$$\begin{eqnarray}&&{\dot{\gamma }}_{{rz}}={\dot{\gamma }}_{{zr}}=\displaystyle \frac{{\rm{d}}{u}_{z}}{{\rm{d}}r},\end{eqnarray} \tag{ 9 }$$

and hence the only non-constant components of stress in all of our models are σ_rz = σ_zr.

These observations, along with denoting the constant pressure gradient, $\partial p/\partial z,$ by −G, reduce the Navier–Stokes equations to

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r{\sigma }_{{zr}}\right)=\rho g-G.\end{eqnarray} \tag{ 10 }$$

Integrating equation (10) with the boundary condition from symmetry that σ_zr(0) = 0, and introducing the wall stress constant,

$$\begin{eqnarray}&&{\sigma }_{{\rm{w}}}=\displaystyle \frac{1}{2}a(\rho g-G),\end{eqnarray} \tag{ 11 }$$

we see that the momentum equation is reduced to

$$\begin{eqnarray}&&{\sigma }_{{zr}}={\sigma }_{{\rm{w}}}\displaystyle \frac{r}{a},\end{eqnarray} \tag{ 12 }$$

where the relation between σ_zr and u_z will depend on which model we choose.

In all our models, velocity boundary conditions on the pipe will be the no-slip condition, and a symmetry observation,

$$\begin{eqnarray}&&{u}_{z}(a)=0,\qquad \displaystyle \frac{{\rm{d}}{u}_{z}}{{\rm{d}}r}(0)=0.\end{eqnarray} \tag{ 13 }$$

3.2.1. Newtonian model

In the Newtonian model, the stress term is

$$\begin{eqnarray}&&{\sigma }_{{zr}}={\mu }_{{\rm{n}}}{\dot{\gamma }}_{{zr}}={\mu }_{{\rm{n}}}\displaystyle \frac{{\rm{d}}{u}_{z}}{{\rm{d}}r},\end{eqnarray} \tag{ 14 }$$

which when substituted into equation (12) and integrated with boundary conditions, equation (13), gives us the velocity profile we expect for Poiseuille flow,

$$\begin{eqnarray}&&{u}_{z}=({r}^{2}-{a}^{2})\displaystyle \frac{(\rho g-G)}{4{\mu }_{{\rm{n}}}}.\end{eqnarray} \tag{ 15 }$$

Note that we should not be worried that r² − a² is negative: in the absence of a pressure gradient G, we would expect fluid to flow down the tube, in the negative z-direction.

With constant values from table 1, figure 3(a) shows this parabolic velocity profile.

Zoom In Zoom Out Reset image size

3.2.2. Power-law model

In the power-law model, the stress term is given by

$$\begin{eqnarray}&&{\sigma }_{{zr}}={\mu }_{{\rm{p}}}{\dot{\gamma }}^{n-1}{\dot{\gamma }}_{{zr}}={\mu }_{{\rm{p}}}{\left(\displaystyle \frac{{\rm{d}}{u}_{z}}{{\rm{d}}r}\right)}^{n},\end{eqnarray} \tag{ 16 }$$

since ${\dot{\gamma }}_{{zr}}$ and ${\dot{\gamma }}_{{rz}}$ are the only non-zero components of the shear rate tensor. When substituted into equation (12), rearranged and integrated with boundary conditions, equation (13), we have the velocity profile

$$\begin{eqnarray}&&{u}_{z}={\left(\displaystyle \frac{{\sigma }_{{\rm{w}}}}{{\mu }_{{\rm{p}}}a}\right)}^{1/n}\displaystyle \frac{{r}^{(n+1)/n}-{a}^{(n+1)/n}}{(n+1)/n}.\end{eqnarray} \tag{ 17 }$$

In the case ${\mu }_{{\rm{p}}}={\mu }_{{\rm{n}}},$ n = 1, this reduces to our Newtonian result. With constant values from table 1, the velocity profile is also shown in figure 3(a).

3.2.3. Casson's model

Substituting equation (12) into Casson's model, we have

$$\begin{eqnarray}&&\sqrt{{\sigma }_{{\rm{w}}}\displaystyle \frac{r}{a}}=\sqrt{{\mu }_{{\rm{c}}}{\dot{\gamma }}_{{zr}}}+\sqrt{{\sigma }_{{\rm{y}}}}=\sqrt{{\mu }_{{\rm{c}}}\displaystyle \frac{{\rm{d}}{u}_{z}}{{\rm{d}}r}}+\sqrt{{\sigma }_{{\rm{y}}}},\end{eqnarray} \tag{ 18 }$$

for σ_zr > σ_y. We solve for u_z by rearranging and integrating (with the boundary conditions, equation (13), again), and, specifying the yield stress, ${\sigma }_{{\rm{w}}}={\sigma }_{{\rm{y}}}a/{r}_{{\rm{c}}},$ we find

$$\begin{eqnarray}{u}_{z}(r)=\left\{\begin{array}{ll}\displaystyle \frac{{\sigma }_{{\rm{w}}}}{2a{\mu }_{{\rm{c}}}}\left(\displaystyle \frac{8}{3}\left({a}^{3/2}-{r}^{3/2}\right){r}_{{\rm{c}}}^{1/2}-2(a-r){r}_{{\rm{c}}}-\left({a}^{2}-{r}^{2}\right)\right) & {\rm{if}}\;{\rm{}}{r}_{{\rm{c}}}\leqslant r\leqslant a\\ {u}_{z}({r}_{{\rm{c}}}) & {\rm{if}}\;{\rm{}}0\leqslant r\lt {r}_{{\rm{c}}}\end{array}\right..\end{eqnarray} \tag{ 19 }$$

Rearranging equation (12) and using values from table 1 gives us an estimate for the critical radius r_c of 2.4 mm, only 12% of the radius of the pipe. This velocity profile has been illustrated in figure 3(a).

The Newtonian model exhibits the well-known parabolic velocity profile. In our shear-thinning models—the power-law and Casson models—the profile is more blunted and eventually a 'plug' of high-viscosity fluid forms in the centre of the channel. The plug flows with nearly uniform velocity, while the velocity gradients are confined to regions close to the wall.

Scaling all the profiles on their maxima to compare their shapes (see figure 3(b)), we see that with small r_c, the Casson profile looks very similar to the parabola seen in the Newtonian case, with a slight flatness at the centre: indeed, equation (19) reduces to the Newtonian expression, equation (15), in the limit ${r}_{{\rm{c}}}\to 0.$ With larger r_c, it resembles much more closely the power-law profile, with the plug flow in the centre. However, in the Casson model, the flatness seen is a genuine plug: it is truly flat, with $\partial {u}_{z}/\partial r=0$ in this area. In the power-law profile, the profile towards the centre is only approximately flat, forming instead a 'pseudo-plug', the flatness of which increases with the strength of shear-thinning of the fluid (as $n\to 0$).

The maximum speed of the power-law fluid is 40% of the Newtonian maximum, which is itself 24% of the Casson fluid. These differences are to be expected and are not features of the models themselves, but rather because the parameters have been sourced from different industrial experiments. These results in the pipe will be compared with those from dome flow in the next section, and the similarity will allow us to disregard Casson's model from here on in, since its profile shape falls between the other two models.

Since we are modelling the dome as the top half of a sphere, we set up spherical coordinates (r, θ, ϕ) where θ is the inclination angle, 0 ≤ θ < π, and ϕ is the azimuthal angle, 0 ≤ ϕ < 2π. We label the radius of the sphere R, and the film thickness h. We expect axisymmetric flow so we can observe h = h(θ).

Again we are attempting to solve the Navier–Stokes equations, equation (1), expecting steady, axisymmetric flow, so we can set u_ϕ, and derivatives with respect to ϕ and time, to zero. Even with the simplifications afforded by this (e.g. ${\dot{\gamma }}_{r\phi }={\dot{\gamma }}_{\theta \phi }=0$) and our constitutive equations (e.g. hence ${\sigma }_{r\phi }={\sigma }_{\theta \phi }=0$), scaling arguments need to be made to attempt an analytical solution.

Scaling is an important technique that gives access to solving difficult problems by providing a rigorous method of deciding when variables are insignificant compared to others. The method is, for each variable, to divide through by a typical (characteristic) value, so that the resultant 'scaled', dimensionless, variables are all of the same order. Then we can compare terms by looking at the sizes of the resulting dimensionless groups of these characteristic values.

We introduce scaled variables (with hats) as follows:

$$\begin{eqnarray}&&\widehat{{u}_{\theta }}=\displaystyle \frac{{u}_{\theta }}{U},\qquad \widehat{{u}_{r}}=\displaystyle \frac{{u}_{r}}{V},\qquad \widehat{h}=\displaystyle \frac{h}{H},\qquad \widehat{y}=\displaystyle \frac{y}{H}=\displaystyle \frac{r-R}{H},\end{eqnarray} \tag{ 20 }$$

where y is a shifted (but still dimensional) radial coordinate which emphasizes the physical nature of the problem: the fluid occupies the space 0 ≤ y ≤ h. U and V are characteristic speeds in the θ- and r-direction respectively, H is a characteristic film thickness, and R is our already-established dome radius.

The continuity equation, equation (8), written in spherical coordinates, is

$$\begin{eqnarray}&&0=\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{u}_{r}\right)+\displaystyle \frac{1}{\mathrm{sin}\theta }\displaystyle \frac{\partial }{\partial \theta }\left({u}_{\theta }\mathrm{sin}\theta \right).\end{eqnarray} \tag{ 21 }$$

Substituting in the scaled variables gives us

$$\begin{eqnarray}&&0=\displaystyle \frac{1}{H\widehat{r}+R}\displaystyle \frac{1}{H}\displaystyle \frac{\partial }{\partial \widehat{r}}\left({\left(R+H\widehat{r}\right)}^{2}\;V\widehat{{u}_{r}}\right)+\displaystyle \frac{1}{\mathrm{sin}\theta }\displaystyle \frac{\partial }{\partial \theta }\left(U\widehat{{u}_{\theta }}\mathrm{sin}\theta \right).\end{eqnarray} \tag{ 22 }$$

Then, since the two terms in the sum balance each other, the characteristic values should balance each other as well. Noting that $H\widehat{r}+R\sim R$ (since $H\ll R$), the characteristic values of both terms, respectively, are

$$\begin{eqnarray}&&\displaystyle \frac{1}{R}\displaystyle \frac{1}{H}{R}^{2}\;V\quad \sim \quad U.\end{eqnarray} \tag{ 23 }$$

Rearranged, this tells us

$$\begin{eqnarray}&&V\sim U\displaystyle \frac{H}{R}\ll U,\end{eqnarray} \tag{ 24 }$$

(since H/R is small) which tells us how these characteristic values interact with each other. In particular, it tells us what we already know: flow is predominantly in the θ-direction.

We then use the same scaling technique on the Navier–Stokes (momentum) equations in spherical coordinates. For example, the θ-momentum equation is

$$\begin{eqnarray}\rho \left[{u}_{r}\displaystyle \frac{\partial {u}_{\theta }}{\partial r}+\displaystyle \frac{{u}_{\theta }}{r}\displaystyle \frac{\partial {u}_{\theta }}{\partial \theta }\right] & = & -\displaystyle \frac{1}{r}\displaystyle \frac{\partial p}{\partial \theta }-\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{\sigma }_{r\theta }\right)+\displaystyle \frac{1}{r\mathrm{sin}\theta }\displaystyle \frac{\partial }{\partial \theta }\left({\sigma }_{\theta \theta }\mathrm{sin}\theta \right)\right.\\ & & \left.+\displaystyle \frac{{\sigma }_{r\theta }}{r}-\displaystyle \frac{{\sigma }_{\phi \phi }\mathrm{cot}\theta }{r}\right]+\rho g\mathrm{cos}\theta .\end{eqnarray} \tag{ 25 }$$

Introducing the scaled variables, we find that the sizes of these terms are

$$\begin{eqnarray}&&\rho \left[\displaystyle \frac{{U}^{2}}{R}+\displaystyle \frac{{U}^{2}}{R}\right]=(?)-\mu \left[\displaystyle \frac{U}{{H}^{2}}+\displaystyle \frac{U}{{R}^{2}}+\displaystyle \frac{U}{{RH}}-\displaystyle \frac{U}{{R}^{2}}\right]+(?),\end{eqnarray} \tag{ 26 }$$

where we have marked by (?) the size of terms that we do not know. Dividing through by μ U/H², the relative sizes are

$$\begin{eqnarray}&&{\varepsilon }^{2}{Re}\left[1+1\right]=(?)-\left[1+{\varepsilon }^{2}+\varepsilon -{\varepsilon }^{2}\right]+(?),\end{eqnarray} \tag{ 27 }$$

where $\varepsilon =H/R\ll 1,$ and ${Re}$ is the Reynolds number: a non-dimensional number which is large in fast, low-viscosity flows; and is small in slow, high-viscosity flows. For a power-law fluid, the Reynolds number, is defined to be

$$\begin{eqnarray}&&{Re}=\displaystyle \frac{\rho {UR}}{\mu }=\displaystyle \frac{\rho {U}^{2}R}{H{\mu }_{{\rm{p}}}}{\left(\displaystyle \frac{H}{U}\right)}^{n},\end{eqnarray} \tag{ 28 }$$

since the viscosity is given by $\mu ={\mu }_{{\rm{p}}}{\dot{\gamma }}^{n-1}.$ Using the values from table 1, we find that for the Newtonian data, ${Re}=0.64,$ and for the power-law data, ${Re}=2.9.$ Both of these are of the order we expect.

In fact, here we have the 'reduced Reynolds number', ${\varepsilon }^{2}{Re},$ which is familiar in thin-film problems. Since the Reynolds number is of order 1, we are able to not only discount those viscous terms on the right-hand side of equation (27) which are multiplied by , but also the inertial terms on the left-hand side.

Using this process on the r-equation at this point, we can deduce the scaling for the pressure term using the same method as with the continuity equation, finding it to be negligible.

Putting this into the θ-equation, we continue our scaling analysis until we are left with

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\sigma }_{r\theta }}{\partial r}=-\rho g\mathrm{sin}\theta ,\end{eqnarray} \tag{ 29 }$$

i.e., the viscosity term balances with the gravity term. This is what we would expect from a fluid with a low Reynolds number.

Recalling ${\sigma }_{r\theta }={\mu }_{{\rm{p}}}{\dot{\gamma }}^{n-1}{\dot{\gamma }}_{r\theta },$ we can substitute this into equation (29) and perform similar scaling arguments to arrive at

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{u}_{\theta }}{\partial {y}^{2}}+\displaystyle \frac{\rho g\mathrm{sin}\theta }{{\mu }_{{\rm{p}}}n{2}^{(1-n)/2}}{\left(\displaystyle \frac{\partial {u}_{\theta }}{\partial y}\right)}^{1-n}=0,\end{eqnarray} \tag{ 30 }$$

with boundary conditions

$$\begin{eqnarray}&&{u}_{\theta }(0)=0,\qquad \displaystyle \frac{\partial {u}_{\theta }}{\partial y}(h)=0.\end{eqnarray} \tag{ 31 }$$

By writing ${U}_{\theta }=\partial {u}_{\theta }/\partial y$ and defining $\kappa =\rho g\mathrm{sin}(\theta )/({\mu }_{{\rm{p}}}{2}^{(1-n)/2}),$ the differential equation becomes separable⁵

$$\begin{eqnarray}&&\displaystyle \frac{\partial {U}_{\theta }}{\partial y}+\displaystyle \frac{\kappa }{n}{U}_{\theta }^{1-n}=0,\end{eqnarray} \tag{ 32 }$$

with solution

$$\begin{eqnarray}&&{U}_{\theta }(y)={[\kappa (h-y)]}^{1/n}.\end{eqnarray} \tag{ 33 }$$

Integrating this with our boundary conditions, (31), gives us the expression for flow down the dome as

$$\begin{eqnarray}&&{u}_{\theta }(y)=\displaystyle \frac{n}{1+n}{\left(\displaystyle \frac{\rho g\mathrm{sin}\theta }{{\mu }_{{\rm{p}}}{2}^{(1-n)/2}}\right)}^{1/n}\left[{h}^{(n+1)/n}-{(h-y)}^{(n+1)/n}\right].\end{eqnarray} \tag{ 34 }$$

Although we have found the expression for the velocity, equation (34), as desired, it relies on knowing the film thickness, h(θ), which is itself a function of angle down the dome. Fortunately we can introduce a constant, measurable quantity which allows to us find the thickness function h(θ). Once we have found this, we can put it back into our velocity expression.

This constant quantity is the flux. For this we have to assume that the flow is steady, i.e. there are no variations in the amount of chocolate being pumped out of the top of the pipe. Then at any time, since the chocolate falls down smoothly, the amount of chocolate passing through all points on the dome with some angle, is the same as the amount passing through any other set of points with another angle. In particular, at any given time, the amount of chocolate being pumped onto the top of the dome is the same as the amount of chocolate falling off the dome.

Mathematically, the flux, Q, is described by the velocity of the chocolate, integrated over the area it is passing through:

$$\begin{eqnarray}&&Q=2\pi \mathrm{sin}\theta {\displaystyle \int }_{0}^{h(\theta )}(R+y){u}_{\theta }(y)\;{\rm{d}}y,\end{eqnarray} \tag{ 35 }$$

where the sin θ term before the integral comes from the 'flat radius' of the dome at any given angle θ, which varies through the chocolate from $R\mathrm{sin}\theta$ to $(R+h)\mathrm{sin}\theta .$

This integration appears unwieldy but requires simply an integration by parts. However, doing this would leave us unable to find an explicit expression for the film thickness h(θ). Instead we once again make the observation that the film thickness is much smaller than the radius of the dome, y ≤ h ≪ R, and so our flux expression becomes

$$\begin{eqnarray}&&Q=2\pi \mathrm{sin}\theta {\displaystyle \int }_{0}^{h(\theta )}{{Ru}}_{\theta }(y)\;{\rm{d}}y,\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&=\;\tau {h}^{(2n+1)/n}(\theta ){\mathrm{sin}}^{(n+1)/n}\theta ,\end{eqnarray} \tag{ 37 }$$

where

$$\begin{eqnarray}&&\tau =\displaystyle \frac{{2}^{(3n-1)/2n}n\pi R}{2n+1}{\left(\displaystyle \frac{\rho g}{{\mu }_{{\rm{p}}}}\right)}^{1/n}.\end{eqnarray} \tag{ 38 }$$

Since we can measure the film thickness at any angle θ with simply a ruler, we choose to measure it at the base of the dome (θ = π/2). Using the measured value of 2 mm, as well as the other relevant parameters from table 1, the Newtonian model estimates a flux of 2.1 ml s⁻¹, and the power-law model estimates 0.015 ml s⁻¹: a considerable difference. Some intuition suggests this these are both too small: we can directly measure the flux by capturing all the chocolate that falls in a 10 s timeframe, measuring the volume, and dividing by 10 s. Doing so gives us a measured flux of 14 ml s⁻¹.

Rearranging equation (37) to find an equation for the film thickness h(θ), having set our constant value of Q to be when θ = π/2, we find that

$$\begin{eqnarray}&&h(\theta )=\displaystyle \frac{h(\pi /2)}{{\mathrm{sin}}^{(n+1)/(2n+1)}\theta }.\end{eqnarray} \tag{ 39 }$$

The simplicity of this result is surprising and satisfying: the film thickness for a power-law fluid at any angle is given by a single trigonometric function, scaled on the (easily measurable) final film thickness. This equation for both models has been plotted in figure 4, with final thickness h(π/2) = 2 mm.

Zoom In Zoom Out Reset image size

We also now use these expressions for film thickness, h, in our expression for the velocity of each fluid, u_θ(y, θ), equation (34). Since the expression for velocity includes rheological parameters (from table 1), in order to compare the two models we have set the flux to be, in both cases, that which we measured in experiment.

The velocity has been plotted as a function of angle, at some film thicknesses, in figure 5; and at a fixed angle, at varying thicknesses, in figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have found results which, from looking at figures 4 and 5 seem sensible. But the severe underestimate of the flux in equation (37) is troubling. Such a low estimate for the flux suggests that the parameters and boundary conditions we have used lead to very slow flow. Indeed, equation (34) suggests a velocity at the upper surface of the chocolate, when it is about to fall off, of 17 μm s⁻¹ for the power-law case, many orders of magnitude smaller than the observed speed of 10 cm s⁻¹.

We should note, though, the sensitivity to our measurements. Our expression for the flux, equation (37), shows that $Q\propto {h}^{(2n+1)/n}.$ In our power-law model, n = 1/3, so we have Q ∝ h⁵: clearly this expression for the flux is highly sensitive to the accuracy of our height measurement. However, the height of the film was measured quite crudely, by placing some paper in the film and measuring the height of the chocolate left on the paper. A measured final height of 1 mm leads to a power-law final velocity of 1 μm s⁻¹, whereas a final height of 2 mm leads to a velocity of 17 μm s⁻¹: a significant difference which we need to be careful about when we have such primitive measurements. Indeed, this is a typical problem when data-fitting for strongly shear-thinning fluids.

A likely inaccuracy is the 'no-slip condition', the first statement in equation (31): it may be that the chocolate is slipping over the dome surface. If we instead knew this speed on the boundary, our theory would provide higher speeds, perhaps better matching the observed flux.

Another source of inaccuracy is our rheological parameters. Despite extracting data from sources for milk chocolate at 40 °C , there is still a large variety of different chocolates, made in different ways and with different proportions of ingredients. The large difference in the two pieces of data in figure 2 is typical of what we find when trying to source rheological data. That the prediction for the fluid is so small, suggests that the chocolate is not as viscous as our data suggests. Indeed, in section 3, figure 3(b) showed that the parameters in the Casson model—a model since discounted as displaying intermediate behaviour between the Newtonian and power-law models—predict much higher speeds in the pipe, which would transfer to higher flow speeds here as well.

Rather than take the data to predict the speed and thickness, it would be nice to be able to measure these to infer the rheological parameters. With accurate measuring equipment, it should be possible to solve equation (37) for the power-law parameters, μ_p and n, if we can measure accurately the height of the fluid in two places on the dome. However, without such equipment to test it ourselves, this is the best we can do, and so it serves as a useful lesson on the uncertainty of rheological measurements.

A final source of error could be that we assume that the chocolate runs steadily down the dome. On closer observation, it tends to come down in waves, as it is pushed out of the top of the pipe. However, this should not account for the wildly small velocity prediction.

Problems aside, let us have a look at the results. The velocity profiles for each fluid from equation (34), when plotted at different film thicknesses in figure 6, match the velocity profiles for half a pipe flow from the previous section. The plug flow seen for the power-law model is also seen in the dome flow. This is not surprising, as when we arrive at equation (29), we see that the curvature of the problem has disappeared. The balance of viscosity and gravity is what we would find for a fluid travelling down an flat inclined plane at an angle, θ, to the horizontal. In other words, to leading order in the film aspect ratio, the fluid does not experience the curvature of the substrate. Although we have modelled the dome as a hemisphere, we could choose any smooth shape for the dome and the velocity profiles would still describe the flow for each locally sloped piece of surface.

The velocity profiles in figure 5 show that the shear-thinning power-law fluid is slower than the Newtonian fluid at varying film thicknesses throughout the majority of flow (matching our observation in the pipe) although close to the surface of the dome, it begins to match it.

In figure 4, we can see that both models give qualitatively the same thickness profile, well up to high enough on the dome where the flow would be affected by the previous part of the fountain. Given the complication of the power-law model for no further insight, use of the Newtonian model for further work here would be sufficient.

The water bell problem was first investigated by Taylor and Howarth [23], where they derived an equation of motion for a water bell trajectory. Assuming the flow to be axisymmetric, we can take a vertical slice of the curtain (see figure 8) and consider the forces acting on a fluid element: gravity, surface tension, and internal and external pressure difference.

Zoom In Zoom Out Reset image size

The form of the gravitational force is as we would expect (mass of the fluid element multiplied by acceleration due to gravity); and the pressure force is given by multiplying the volume of the element by the pressure difference, set in the direction normal to the curtain. The force due to surface tension is given by the Young–Laplace equation, again multiplied by the volume of the element and set in the normal direction.

By balancing these forces and resolving (with suitable non-dimensionalization in the style of the previous chapter), Taylor and Howarth found a relationship between the velocity, u(z), and distance of the sheet from the pipe in the centre of the bell, r(z), both functions of vertical distance below the bottom of the dome, z. Some later trigonometric simplifications by Brenner [5, equations (13)–(15)] lead to

$$\begin{eqnarray}&&{r}^{\prime\prime }(u-r)+\left(1+r{^{\prime} }^{2}\right)\left(1+\displaystyle \frac{\beta }{u}r^{\prime} \right)-\alpha r{\left(1+r{^{\prime} }^{2}\right)}^{3/2}=0,\end{eqnarray} \tag{ 40 }$$

the process of which is explained fully in Button [7, ch 5].

This equation has some given initial conditions: r(0) = r₀, the scaled radius of the dome; and r'(0) = v₀, specifying the angle at which the curtain begins its fall. The problem depends on two physical parameters: α accounts for the effect of the inside–outside pressure difference, and β accounts for gravity. They are discussed more in section 5.3.2, but are given by

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{\rho {{Qu}}_{0}{\rm{\Delta }}p}{8\pi {\gamma }_{{\rm{s}}}^{2}}\quad {\rm{and}}\quad \beta =\displaystyle \frac{\rho {gQ}}{4\pi {\gamma }_{{\rm{s}}}{u}_{0}},\end{eqnarray} \tag{ 41 }$$

where we can see they are composed of terms we would expect from the three forces we are balancing: γ_s is surface tension; Δp is the pressure difference between the inside and outside of the sheet; and g is the acceleration due to gravity. We also have u₀, our initial velocity; and Q, the (constant) flux.

We assume that the sheet is sufficiently thin so that we can assume the velocity to be constant through a cross-section at a given height. Since the surface of the liquid film is itself a streamline and this model is for inviscid fluid, we can use Bernoulli's equation for any arbitrary point along a streamline where gravity is constant:

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\rho {u}^{2}-\rho {gz}+p={\rm{constant}},\end{eqnarray} \tag{ 42 }$$

noting that we have switched the sign on the $\rho {gz}$ term from the usual convention, as we are measuring z downwards.

Since the streamline is also a free surface, the pressure, p, on it is constant, so it can be absorbed into the right-hand side. With the non-dimensionalized boundary conditions u = 1 at z = 0, Bernoulli's equation reduces to

$$\begin{eqnarray}&&{u}^{2}=1+2\beta z,\end{eqnarray} \tag{ 43 }$$

which expresses the increase of momentum in the direction of the flow due to the acceleration of gravity.

We are left with the complete governing equations, equations (40) and (43), to solve. Analytically this is only possible where α = β = 0. Since we do not expect any significant pressure difference between the inside and outside of the sheet, we can set α = 0. However, we cannot say this for β, since this is the Bond number (or Eötvös number), which represents the balance of gravitational forces to surface tension forces which we are left with from figure 8. With values from table 1, we find β = 6.3.

An approach we can take is to consider the β ≫ 1 limit, where the the flow is dominated by gravity rather than surface tension. In this limit, the curtain, starting at r(0) = r₀ (the scaled radius of the dome), falls further downwards than it does inwards, and we can express r as a regular series in $\delta =1/\beta \ll 1$

$$\begin{eqnarray}&&r(z)={r}_{0}+\delta {r}_{1}(z)+{\delta }^{2}{r}_{2}(z)+\cdots .\end{eqnarray} \tag{ 44 }$$

Differentiating, we have

$$\begin{eqnarray}&&{r}^{\prime }(z)=\delta r{{}^{\prime }}_{1}(z)+{\delta }^{2}r{{}^{\prime }}_{2}(z)+\cdots \end{eqnarray} \tag{ 45 }$$

and substituting this into equation (40) and taking the leading order terms, we have

$$\begin{eqnarray}&&1+\displaystyle \frac{r{{}^{\prime }}_{1}(z)}{u(z)}=0.\end{eqnarray} \tag{ 46 }$$

Rearranging for $r{{}^{\prime }}_{1}(z)$ and taking the leading order term from equation (45), we are left with

$$\begin{eqnarray}&&{r}^{\prime }(z)=-\displaystyle \frac{u(z)}{\beta },\end{eqnarray} \tag{ 47 }$$

which we can substitute into equation (43).

Integrating the resulting equation with the boundary condition r(0) = r₀, we have an expression for the sheet profile r(z),

$$\begin{eqnarray}&&r(z)={r}_{0}+\displaystyle \frac{1-{(1+2\beta z)}^{3/2}}{3{\beta }^{2}}.\end{eqnarray} \tag{ 48 }$$

This equation (in dimensional form) has been plotted in figure 9, with data from table 1.

Zoom In Zoom Out Reset image size

Finding an equivalent profile for even a Newtonian viscous liquid requires numerical work, and we were pleased to find Gilio et al [12], where they study both inviscid and viscous two-dimensional, falling, nonsteady liquid sheets. They derive governing equations using a different method, and solve their viscous equations numerically, using finite-difference methods.

In their viscous model they use a fluid with Newtonian viscosity 10 Pa s, close to our chocolate value of 14 Pa s. We were able to manipulate our inviscid graph, from equations (40) and (43) to fit theirs, to find an equation governing the viscous velocity.

We do this by scaling the velocity plots from [12, figure 6] to match our inviscid results in equation (43). The viscous profile is remarkably linear and can be extracted from their graph, in non-dimensional form, as

$$\begin{eqnarray}&&u(z)=1+2.12z.\end{eqnarray} \tag{ 49 }$$

We can substitute this velocity into equation (47) to find an equation governing the position of the viscous sheet, namely

$$\begin{eqnarray}&&r(z)={r}_{0}-\displaystyle \frac{z+1.06{z}^{2}}{\beta }.\end{eqnarray} \tag{ 50 }$$

In this way we now have a viscous prediction for the profile of our falling sheet, and this has also been plotted (in dimensional form) on figure 9.

The inviscid velocity profile, equation (43), is somewhat validated by the results of Brunet et al [6]. They use a silicon oil with similar density and surface tension to chocolate (ρ = 970 kg m⁻³, γ_s = 20.4 mN m⁻¹), but with a viscosity of 0.2 Pa s, which is smaller than that of chocolate (but considerably higher than water). They find by direct measurements that equation (43) is a reasonable approximation of their results.

Comparing the inviscid sheet profile with the observed falling profile in figure 9, we see that the sheet falls inwards due to surface tension with slope of the same order we expect. However, it predicts that the sheet falls in about 3 cm over the 7 cm drop, when we measure in experiment about half of that. Secondly, in our observations, the sheet starts falling at a slant, rather than downwards, and we have lost the ability to set r'(0) in performing the perturbation analysis in equation (5.1), which reduced the order of the governing equation. The viscous sheet prediction, however, falls 1.5 cm inwards, better matching our observations.

5.3.1. Validity of inviscid approximation

Having derived profiles from two models, we now show that we expect the inviscid model to fare poorly, since the viscous effects are expected to be significant. We compare the total energy in the system with the amount used in the work done by viscous forces in the falling fluid. We find that even in the least-possible case for the viscous forces, they are of comparable size, meaning that ignoring viscous forces, as we have done in our inviscid approximation, is unlikely to give us good results.

We start by finding the total energy in the system, i.e. the sum of kinetic and gravitational energy. The velocity in our sheet is predominantly vertical, occupying in cylindrical polar coordinates the (dimensional) region V, given by R ≤ r ≤ R + h, 0 ≤ θ < 2π, 0 ≤ z ≤ ℓ, recalling that the z-axis points downwards. Then the kinetic energy E_K in the system is given by

$$\begin{eqnarray}&&{E}_{{\rm{K}}}=\displaystyle \frac{1}{2}\rho {\displaystyle \int }_{V}{u}^{2}\;{\rm{d}}V=\displaystyle \frac{1}{2}\rho A\left({u}_{0}^{2}{\ell }+g{{\ell }}^{2}\right),\end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray}&&A={\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{R}^{R+h}r\;{\rm{d}}r\;{\rm{d}}\theta =\pi \left[{(R+h)}^{2}-{R}^{2}\right],\end{eqnarray} \tag{ 52 }$$

the cross-sectional area of a horizontal slice of the curtain.

The gravitational potential energy E_P of the system is given by

$$\begin{eqnarray}&&{E}_{{\rm{P}}}=A{\displaystyle \int }_{0}^{{\ell }}\rho {gz}\;{\rm{d}}z=\rho {gX}\displaystyle \frac{{{\ell }}^{2}}{2}.\end{eqnarray} \tag{ 53 }$$

Therefore the total energy instantaneously in the falling sheet, when we plug in the parameters from table 1, is E_K + E_P = 0.044 J.

We want to compare this with the amount of energy dissipated as the sheet falls. Cauchy's energy equation for the rate of energy dissipation, K, due to viscosity in a viscous Newtonian fluid of volume V and viscosity μ is given by [15]

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}K}{{\rm{d}}t}=-{\displaystyle \int }_{V}4\mu {\dot{\gamma }}^{2}\;{\rm{d}}V.\end{eqnarray} \tag{ 54 }$$

Here, $\dot{\gamma }$ is the scalar shear rate, first introduced in equation (3). Given that the velocity in our sheet is predominantly vertical, we say that the dominant rate of strain tensor component is ${\dot{\gamma }}_{33},$ given by du/dz. In dimensional form, equation (43) is given by

$$\begin{eqnarray}&&{u}^{2}={u}_{0}^{2}+2{gz},\end{eqnarray} \tag{ 55 }$$

which we can differentiate to find, substituting in ${\dot{\gamma }}_{33}$

$$\begin{eqnarray}&&{\dot{\gamma }}_{33}=\displaystyle \frac{g}{{({u}_{0}^{2}+2{gz})}^{1/2}}.\end{eqnarray} \tag{ 56 }$$

Incompressibility in the fluid (${\boldsymbol{\nabla }}\cdot {\bf{u}}=0$) requires that the shear rate tensor is traceless. We then have some flexibility with values of the other diagonal rate of strain elements ${\dot{\gamma }}_{11}$ and ${\dot{\gamma }}_{22}.$ The least energy change rate in equation (54) occurs when ${\dot{\gamma }}_{11}={\dot{\gamma }}_{22}=-\displaystyle \frac{1}{2}{\dot{\gamma }}_{33},$ and the most when ${\dot{\gamma }}_{11}=-{\dot{\gamma }}_{33}$ and ${\dot{\gamma }}_{22}=0.$ Given that we are looking for a lower bound on the energy dissipation, we take the least-case scenario and after substituting into equation (54), we get

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}K}{{\rm{d}}t}=-{\displaystyle \int }_{V}3\mu \displaystyle \frac{{g}^{2}}{{u}_{0}^{2}+2{gz}}\;{\rm{d}}V=-\displaystyle \frac{3}{2}\mu {gX}\mathrm{log}\left(\displaystyle \frac{2g{\ell }}{{u}_{0}^{2}}+1\right).\end{eqnarray} \tag{ 57 }$$

Again with the parameters from table 1, we find that in a tenth of a second, an estimate for how long it takes chocolate to complete its fall, in the least-case scenario, the energy dissipated by viscous forces K is 0.046 J.

We have shown that, if we use an inviscid model to predict the fluid velocity of our viscous fluid, the dissipation we have neglected is as large as the total energy in the system. Though our energy estimates are rough, it is clear that the neglect of viscous effects is not justifiable, and we can expect the inviscid approximation to give poor predictions.

5.3.2. Importance of physical parameters

The governing equation for our sheet position and velocity profiles, equation (40), relies on the physical parameters α and β, given in equation (41).

Since the pressure difference is negligible, given that the falling sheet never forms an airtight seal, we set α = 0, but we can see that the Bond number, β, depends on physical parameters of the chocolate.

Figure 10 shows the profile of the falling sheets in the inviscid and viscous cases for substances with different values of β. Lower values of β, corresponding primarily to increased surface tension since the densities are of similar order, bring the sheet in further. Interestingly, it affects the viscous sheet more than the inviscid sheet, with the viscous sheet falling further inwards in the higher surface tension case than the inviscid prediction.

Zoom In Zoom Out Reset image size

From our analyses in the three separate domains of the fountain, we have found how the chocolate moves round the fountain under different models, and were able to compare them to what we see with our actual fountain. In short:

• Pumped pipe flow. For power-law flows, plug flow is expected in a pipe with a pressure gradient. For Newtonian flows, a parabolic flow is expected. The analysis here allows us to narrow our choice of models to just two.
• Dome flow. Gravitational and viscous effects balance, and the dome shape is not important beyond the local slope. The chocolate thins and slows down as it travels over the dome. The Newtonian model performs as well as the power-law model in predicting the thickness of the chocolate.
• Curtain flow. Surface tension is the dominating factor in pulling the sheet inwards. An inviscid model is not expected to predict well how far inwards the sheet will fall, but a viscous one would.

We have taken a novelty, engaging consumer item and presented a project which uses it to introduce the world of non-Newtonian fluid mechanics. Moreover, we have found that it works particularly well in bringing this exciting area of research to many different levels of technical ability.

The process of deriving the key results is necessarily mathematical, and highly suited to a final-year undergraduate course or project. Breaking the fountain into its three natural components, we see that each serves a distinct and valuable educational purpose.

• Pumped pipe flow. Pressure-driven flow in a pipe is a simple flow problem, ideal for developing the theory of non-Newtonian fluid flow and investigating the phenomenology of different model fluids. Students will have met this problem for Newtonian flow before in a fluids course, so it serves as a reminder as well as an interesting extension.
• Dome flow. Thin-film flow over the dome provides an introduction to lubrication theory, by the use of scaling, an important technique which may have been seen already in a different area of applied mathematics. With very mild constitutive assumptions, analytical solutions can be calculated for power-law as well as Newtonian fluids.
• Curtain flow. Finally, the curtain or sheet region in which the molten chocolate falls freely under gravity is intractable analytically (and beyond the scope of an undergraduate level project to solve numerically). However, perturbation analysis allowed us to find a first prediction. After this, we moved onto finding and interpreting relevant existing literature, working out how to convert existing numerical results to our situation by scaling for our values of viscosity, film thickness and speed. At a research level, the falling sheet leaves the greatest opportunity for further work.

Although the level of mathematics required to derive this is high, the results, however, we have found to be very much of interest to the general public in our engagement work. That general foodstuffs can be divided into different rheological classes based on whether they are shear-thinning, shear-thickening, or simply Newtonian, is interesting; but that we can describe them all mathematically with the power-law model, simply by changing the parameter n, shows the power of mathematics. We have found that showing audiences the effect of varying n on viscosity–shear rate graphs allows them to appreciate why mathematics is important in describing the world. Also, that surface tension brings the chocolate fountain inwards as it falls, settles a debate which chocolate enthusiasts have long wondered about.

The subject matter is intrinsically attractive (tasty, even!) and we hope that others get the opportunity to learn about this area of mathematics with their own chocolate fountain.