A dynamic in-silico model captures the kinetics of 1-d gravity driven instabilities, in gravity or centrifuge, of fluid-infiltrated poroelastic media in a partial differential equation (pde). The pde yields the porosity profile over height and time for the given initial and boundary conditions, during slow isotropic isothermal compaction in counter-current fluid drainage.
Processes included are sedimentation, creaming, caking, synaeresis, (ultra) filtration or settling of sludge, pulp, and slurry, washing, sand and ice cream sculpting, filter pressing, injection moulding, consolidation and subsidence in civil engineering, in gas and oil recovery and in carbon dioxide sequestration, dynamics of biological tissue, etc. The dispersed medium can be anything from (even partially aerated) porous, fibrous, granular dispersions, suspensions, emulsions, fat crystal networks, clothes, soil, rock, or their mixes, and the fluid can be either oleic or aqueous. Typical products are pourable dressings, spreads, mayonnaise, jam, butter, margarine, ice cream, frozen drink, custard, mustard, cheese, ketchup, pudding, yogurt, rock, sand, soil, clothes, muscle tissue, shampoos, creams, soaps and gels, et cetera.
The most important limiting prerequisite is that the incompressible dispersed medium is sufficiently structured and/or concentrated that it compacts isotropically during drainage, without segregation in sizes or in components. The flow rates are so slow that Darcy’s law applies for flow through the deforming porous media and the Newtonian fluid’s flow is Stokian. Carman-Kozeny’s equation applies for permeability. Van Wyk’s equation applies for elasticity. The pde uses four dimensionless Pi numbers. Two Pi numbers are proportional scaling of time and height, one is the initial porosity, and the last Pi number is the ratio of buoyancy and elastic force.
For Unilever, modeling of gravitational instability of products is important to quantify or extrapolate long time behavior during shelf life or use centrifuge data to quickly predict long term shelf performance of products.
For now, the profile is created in a 1-d vessel or reservoir, e.g. vessel or an aquifer, of fixed volume. The medium has initially constant composition and porosity, without gravity or centrifuge effects, as if the medium is freshly mixed or in zero gravity. The boundary conditions (bc) are no flow through the boundaries.
After the compaction process starts, the vertical profile almost immediately converts into three zones separated by two fronts: a clear (depleted) supernatant fluid layer, which borders the initial mix centrally, which borders the caked (elastically or plastically compressed) layer at the other side of the vessel. On which side the fluid layer occurs, depends on the density difference between the fluid and dispersed particulates. In sedimentation, the dispersed solids are heavier, causing the fluid layer to build on top; in creaming, the reverse occurs. The initial dispersive mix in the middle moves at constant speed and constant porosity, subject to buoyancy, and decreases in length as the clear fluid layer and caked layer increase in size. The cleared fluid layer height increases first linearly in time, until the two fronts merge into one front that slows in time due to the elastic force that builds on compaction and counteracts the buoyancy (see Figure 1), until a no-flow balance is reached at equilibrium. As fluid displaces dispersed material in pure counter-current flow, at each level the volume rate of fluid is equal to the volume rate of the dispersed material, but in opposite directions.
Other boundary conditions will create different profiles, for example, during a clothes wash process, in synaeresis, or in filtering compaction, et cetera.
The continuous incompressible fluid (aqueous or oleic) is assumed Newtonian.
Figure 1. Mayonnaises after many weeks on the shelf. On the left an insufficient stabilized product, creamed with a clear layer on the bottom. On the right a better stabilized product.
The isotropic dispersed medium is concentrated or structured in such a way that all incompressible dispersed ingredients move together in the opposite direction of the fluid, while compacting during fluid drainage (counter-current flow). The dispersed medium can be any form, including dispersions, emulsions, clothes, soil, rock, or their mixes, and in any state (aerated, porous, fibrous, granular, et cetera). The fluid wets the dispersed medium. Flow can be driven by buoyancy, centrifuge, applied pressure, simple fall impact of a wetted cloth, et cetera. Temperature is assumed constant. Coarsening (Ostwald ripening of dispersed material and bubbles) is assumed absent. The flow rates are so slow that the fluid follows Stokes flow and Darcy’s law applies for the flow through the deforming porous media. As we are considering the flow in a vessel/reservoir with a closed bottom and closed top (which can be a meniscus), the total sample volume is conserved. The vessel is assumed “wide”, such that any potential side-wall drag is not interfering with the flow rates, and the pore averaged flow is therefore assumed purely 1-dimensional in a vertical upward or downward direction. The total volume below any level does not change: the Darcy speed of particulate downward in sedimentation is therefore at any level exactly equal to the Darcy rate of fluid upward (or both zero). The same applies in reversed mode (e.g. creaming). For the small air bubbles, we assume that the bubble volumes are preserved, which is a simplifying approximation. Volume changes of bubbles are small anyway in the relatively small pressure gradients. The total gas volume fraction must also be constant, as all other ingredients are assumed incompressible and total volume is preserved. Whether the dispersed phase acts elastically or plastically is irrelevant in slow 1-d sustained compression: both for snow and cotton we feel in compression tension the reaction force, but on release of tension snow does not relax back (e.g. is plastically deformed), while cotton expands back (e.g. is elastically deformed).
The fluctuating pore-level Navier-Stokes equations for the fluid flow at sub-pore level are averaged to macroscopic flow using poroelastics. This leads to a Darcy type of fluid flow equation driven by buoyancy forces, slowed down by counteracting elastic forces when the deformable porous medium compacts. The flow equation, together with the volume balance, lead to the partial differential equation (pde), that needs to be solved numerically for the proper initial and boundary conditions. Scaling relations are used to express the pde in scaled dimensionless variables such that one solution actually covers a whole range of actual conditions, dimensions, timings, and variation in parameter values, which cover many processes in various formulations and products.
The pde uses four Pi numbers, based on dimensionless combinations of one value of the seven relevant physical parameters (permeability, elasticity, density, viscosity, height of the vessel, acceleration of free fall, and initial porosity), to characterize any porosity profile as a function of height and time; from the initial waiting and transient to ultimate equilibrium. Two Pi numbers are scaling of time and height, one is the initial porosity, and the last Pi number is the ratio of buoyancy and elastic force.
2. Flow Averaging in Multiphase Flow
The local balances of mass, volume and momentum of each phase, e.g. fluid phase and dispersed structured solid (fabric) phase, are spatially fluctuating functions on the sub-pore level. These equations are averaged to a larger volume, to damp out the local pore-level fluctuations, and find the equations for the mass balance and flow at the pore-averaged level, as described in      . In our case, that is the averaged vertical flow of incompressible fluid in our sedimenting (or creaming) system, described by the local volume balance equation and the local isothermal pore-average momentum balance for slow Darcy volume flow q of an incompressible fluid inside a structured porous medium, driven by buoyancy and elastic forces with , where k is permeability of the porous medium, μ the viscosity of the fluid, Δρ the averaged density difference between the dispersed solids and the fluid, the porosity, and ε the averaged elasticity (compression modulus) of the porous medium under compression. The fluid flow q is the Darcy fluid volume flow per unit area per unit time into the fluid saturated structured dispersed porous medium. If needed, more details regarding derivation and implicit assumptions can be obtained from the author. Combining, we find the pde .
2.1. Darcy Equation
The Darcy equation was derived originally for fluid flow in stagnant soil or rock as . In our case both the fluid and dispersed material flow, hence Darcy’s equation must be adapted. For us, as a geostationary observer, the flows are counter-current with equal flow rates q. If we would be an observer traveling together with the dispersed solids, we would conclude that the material appears stagnant and that therefore Darcy applies exactly. For such an observer the fluid flows with twice the flow rate in apparently stagnant dispersed material, hence Darcy gives . This explains the factor 1/2 in front of Darcy in counter-current flow; e.g. . NB. Both observers do not accelerate relative to each other, thus measure the same forces and pressures.
2.2. Density, Permeability, and Elasticity Characterizations
The fibrous network experiences a buoyancy force per unit volume that is a summation over all non-fluid ingredients moving together in the network with the same speed. The component averaged density difference is then: , which is a constant during the compaction. Here, is the density difference between the density of the particulate i and the density of the fluid (e.g. the buoyancy of the particulate), and is the volume fraction of that particulate, and the summation i is over all the captured materials (structurant, solids, oil droplets, air bubbles, fibres, et cetera) in the fluid. If the non-fluid material is effectively heavier than water, sedimentation occurs, or creaming occurs if lighter. Both permeability k and elasticity ε are functions of the porosity, thus vary over distance and time, and make this second order pde nonlinear and awkward to solve. Cleverer combinations of these two varying parameters can be derived to create better constants and make life easier. The permeability is a function of both porosity and surface area per unit volume. For sphere packings and less regular assemblies of granular media, Carman-Kozeny found a typical experimental correlation , where a is the grain surface area per unit of grain volume and ac an experimental correlation coefficient   . We may assume that the specific surface area a of the dispersed phase does not change significantly during compaction, e.g. during compressional porosity changes. This is a good approximation even for the lower porosities reached in a compacting porous cake, hence we may write , where k0 is the permeability at a reference porosity . The elasticity of a fibrous network pack varies with compression. To describe the elasticity, we will use the Van Wyk’s heuristic equation  . For a fibrous fabric, it relates the applied pressure P to the inverse cube of the volume V, the intrinsic uncompressed volume V00, corrected for the incompressible (high pressure solid limit) volume Vs as follows   : . The equation is independent of the fibre diameter or elasticity, but includes the Young’s modulus E, the mass of the fibres m and the bulk density ρ at low pressure. Further, the equation comprises a dimensionless constant Kc characterising the fibres, of order 0.01 that will vary with fibre orientation. Now, the overall elasticity of the gel network is defined by , thus using the modified Van Wyk’s equation . The fluid-filled pore space is . In vertical compression, the volume of particulates is conserved; . Thus, . So, between two states in compaction . We selected the reference state as the porosity condition, where we can measure the value of porosity dependent physical parameters like k0 and ε0. In practice, this is the initial condition. In the uncompressed state, the pressure is zero when no external forces, like gravity, are applied. In the initial condition, the pressure might not be zero, even without external forces. The material might tend to swell (if excess free fluid is available) or crimp in all directions, even without external force.
If the dispersed material swells initially in an exponentially decaying process, the initial fluid created by the gravitational compaction of the depleted layer is resorbed by the swelling, until the rate of swelling is lower than the compaction rate and the depletion layer starts to build. This leads to an effective extra waiting time before the observable depletion layer (supernatant layer) starts to build. If the dispersed material crimps initially, this can be interpreted, and modeled, as a synaeresis process with an evenly distributed non-zero pressure initially. We will for now assume that sufficient concentration of dispersed material is present initially to avoid crimp.
2.3. The Poroelastic PDE and Its Scaling
For our pde we may write , or using Carman-Kozeny and Van Wyk: , where we introduced the height of the aquifer or vessel H as spatial scaling. If we now introduce a dimensionless scaling parameter , e.g. negative in sedimentation, a dimensionless time , and a dimensionless scaled vertical spatial coordinate with , then pde becomes fully dimensionless as . If we define a dimensionless fluid flow rate as , the dimensionless volume balance becomes , where . We also define a dimensionless height of the supernatant Θ as Θ = Ω/H.
2.4. Scaling Consistent with Buckingham Pi Theorem
The dimensionless pde is consistent with the Buckingham Pi theorem. Next to the dimensionless dependent variable , we have 7 independent parameters [ , Δρ, g, H, k, ε, μ] that are functions of the three base units: length m; weight kg; and time s. Hence, we have 7 − 3 = 4 dimensionless relations between these parameters: in our case the four dimensionless numbers , , and determine the dependent variable . Because is dimensionless, it is a Pi number in itself. As z and H have the same dimension, their ratio is a Pi number  . Thus, we need only one consistent set of the 7 independent parameter values [ , Δρ, g, H, k, ε, μ] to determine the value of the 4 independent dimensionless Pi numbers which characterize the porosity profile over time, from initial waiting, during linear growth and slow down of the clear layer height, to final equilibrium. Alternatively, these parameters like permeability and elasticity, can be determined indirectly from the experimentally characterized profile, for instance, by measuring the height Ω of the clear supernatant layer as a function of time.
The Froude number F is the ratio of the inertia force over gravity force, for a system characterized by a density ρ, a length L, and a speed v: . In the same way, the Cauchy number is the ratio of the inertia force over the elastic force . The F and C numbers are taken squared to linearize in v. Hence, the ratio of the Cauchy number over the Froude number gives the ratio of the buoyancy over the elastic forces: , or in our case more specifically . This is proportional to our . The δ is therefore also a measure of the ratio of gravitational and elastic force. For = 0.3 we expect a balance when δ = 4.2.
2.5. Initial and Boundary Conditions
The initial condition is (for now) constant porosity at t = 0. In practice, any monotonic porosity profile, constant or increasing towards the depletion side, can be taken as the initial start profile. The pertinent boundary conditions (bc) are no flow over the boundaries of the vessel Q = 0 at ξ = 0 and ξ = 1 at all times, hence conservation of dispersed material in the vessel: . Since flow over the boundaries is zero, there can only be a vertical redistribution of the porosity profile internally in the vessel or reservoir. At the cake side boundary, the solid material accumulates, thus the porosity drops, e.g. , which means that in view of the volume conservation that at that cake boundary. If the pde prevails up to the boundary, and can be differentiated twice, then at that no-flow boundary , or . This is certainly so at the cake side boundary. This is a nonlinear Robin type of boundary condition. At the depletion (supernatant) side, the porosity rises in the beginning from towards ® 1, so there is and thus . Thus, as long as the porosity function is not saturated ( < 1), the no-flow boundary condition at the depletion side is and thus . This is the nonlinear Robin type of boundary condition again. In the experiment, the time needed for rising porosity to unity at the depletion side is an effective part of the waiting time: no real visible depletion (supernatant) layer is formed during that period of time.
As long as there is no shock front formed ( ), the porosity obeys the pde over the entire range, hence , e.g. , confirming the Robin nonlinear boundary conditions at both sides.
However, as the porosity reaches = 1 at an infinite thin depletion/supernatant layer, the boundary condition is a no-flow condition at the bottom flank of the depletion layer of thickness Θ, where = 1, while the flank forms a moving boundary (e.g. a type of Stefan problem), a sharp moving shock front, described by Rankine-Hugoniot shock conditions. In the depleted (supernatant) layer is = 1 constantly. The lower boundary of the depleted layer is an almost vertical moving front ending in a sharp nick towards = 1. We may interpret this moving boundary as a space shifting of a Dirichlet boundary condition = 1 (over a spatial range with an extend Θ from the top of the vessel).
If we write the pde as , we recognise the general shape of a nonlinear degenerate elliptic-parabolic pde, the Richards’ equation. The left-hand-side (LHS) is the hyperbolic sedimentation equation for the non-diffusive/non-dispersive sedimentation process. The settling velocity or conservative flux convection coefficient is an increasing function of . The LHS leads therefore to a shock front (kinematic wave) propagating with a speed U( ); the sharp front separates the homogeneous suspension of initial porosity from the clear fluid (the depleted, supernatant layer). Hydrodynamic dispersion modifies the above picture and spreads out the front. The ultimate balance of self-sharpening and dispersion depends on δ, but is difficult to predict a priori, since the pde is nonlinear. The nonlinear nature rules out analytical techniques and consequently the equation can only be solved numerically. See, for instance, the work of Concha and Bustos, and of Bürger, Karlsen, and Evje, on sedimentation and flocculation, around the turn of the century, that describe flocculated suspensions.
One important characteristic is that the effective dispersion coefficient becomes infinitely small for high values of approaching unity. This means that the dispersive broadening becomes very small for the higher porosity values near the start of the cleared zone. We know from the boundary condition that the flank of the cleared zone is also very steep for high values of . This makes the flank sharp and steep, as observed experimentally. On the other end, near the base of the caked zone at ξ = 0, the porosities become small. This means that the caked porosity profile is reasonably flat according to the boundary condition, helped by the fact that the dispersion is large, mixing the porosities over a wide range of ξ (similar to a “random walk”). This flatness of the cake is also observed experimentally.
Basically, sharpening and dispersion are opposing effects on the front, thus they can balance each other, leading to a stabilized dispersed profile that, after a transient time, will propagate with a constant velocity U, without changing its shape  . In the experiments, stable and sharp moving fronts were observed.
Another important characteristic of the porosity profile in gravitation occurs, when starting with a profile that is monotone, neither decreasing (in sediment) nor increasing (in creaming), it will retain that same monotonic characteristic up to equilibrium. Any time or spatial oscillations indicate numerical instabilities in the numerical procedure.
One solution of the pde delivers a full set of curves that cover variations in all parameters for a given initial porosity and for a given ratio δ of buoyancy and elastic force. Repeating the calculations for other values of and δ give all possible solutions that exist in nature for slow counter-current compaction of incompressible structured or concentrated particulates and small gas/air bubbles in an incompressible Newtonian fluid. For each scaled non-dimensional curve (T, ξ) determined by a combination ( , δ), there is in real space and time a fan of curves (t, z) as determined by the proportional scaling of time dimension t = τT and of space dimension z = Hξ. Sedimentation (δ < 0) and creaming δ > 0) are mirror images in vertical direction around the = value when acceleration g is a constant.
2.6. Centrifuge: Radial-Symmetric Flow in Axial Rotation Symmetry
We describe the flow in a centrifuge in time t and in cylindrical coordinates (r, q, z). In a fast spinning centrifuge, the flow is radial symmetric in axial symmetry, hence flow and porosity are a function of time and of radial coordinate r. For flow in radial direction in radial symmetry with . In the gravitational field, the acceleration , is practically constant (g = −9.8125 m/s2), but in a centrifuge the centripetal acceleration is a function of r: , with , where ω the cyclic frequency. The buoyancy force increases pressure with distance from the axis in radial flow for axial symmetry in the centrifuge. In this radial flow in axial symmetry, the volume balance becomes: . For a small rectangular vessel at larger distance from the axis (say 4 times its length), the volume balance equations can be approximated as , where is an average value of the distance of the sample from the rotation axis. We define Z for the sedimentation of a small sample in the centrifuge as . In sedimentation, r0 is the outer radius of the sample and . In creaming, r0 is the inner radius of the sample and and , with Δρ being negative. So, for such a small vessel spinning at large distances from the rotation axis and . So has for the small sample in the centrifuge the same role as ggrav, but often represents a much larger gravity in fast spinning centrifuges: where RGF a large positive number, or in such centrifuges: . If we represent earth’s gravity data as a function of z and compare them with centrifuge data for small samples far from the axis as a function of Z, in for instance a sedimentation experiment, the graphs would be (nearly) identical at the same imposed gravity forces. This also applies when scaled dimensionless, using the same four Pi parameters H, τ, , and δ.
2.7. Equilibrium Profile
In the ultimate equilibrium, the fluid flow over the entire vessel, hence . Integration yields the equilibrium profile , where some reference point. If we take the reference point at the point where in equilibrium: , so at the bottom ξ = 0 where : . The supernatant is a clear fluid layer with = 1 over scaled distance .
In radial symmetry in equilibrium
leads to . At the axis r = 0, the porosity is (virtually) maintained at the initial value as there is still the initial pressure irrespective of rotation speed, and from the axial symmetry, hence const = 0, which leads to the same condition as in constant gravity, but now z is replaced by r: e.g. at equilibrium, using Van Wijck in axisymmetric counter-current flow of incompressible fluid and particulates in a closed vessel in the centrifuge: . The equilibrium profile is over a finite part of the vessel: the rest is depleted zone (in the centrifuge at the inner part of the vessel in sedimentation, at the outer part in creaming), such that the particulate volume initially present over the entire vessel is preserved and compacted over the shorter equilibrium profile in the caked compaction zone. In the strong gravity field of the centrifuge, there is significant compaction, even near the top of the cake, which seems counter intuitive, as there is no push from excess weight above. However, the buoyancy force is present, even at the top of the cake.
2.8. Weak and Strong Form
The strong form requires the function to be continuous and at least twice differentiable over its domain. The solution of our pde: may contain discontinuous jumps, which are weak, around which the function remains continuous, twice differentiable and is therefore strong. If we integrate once over space: . This integral describes the same physics but needs to be only once differentiable. This is a weak formulation and is approximate. . A finite difference formulation is also an approximate and weak solution of the full equation and boundary conditions. When the spacing is smaller, the solution gets stronger. The position of the shock is described by the displacement from the flow from both the upstream and downstream sides, e.g. Rankine-Hugoniot jump conditions, as imposed by the requirement of conservation of energy, momentum, mass and volume.
2.9. Initial Speed of Shock Front
The volume balance can be rewritten as a weak shock front speed , where . Over a shock front the volume balance leads to Rankine-Hugoniot condition . For the initial linear displacement cleared front Θ between and depleted layer with 1 with respectively and , we find . Since the dimensionless clear zone length in sedimentation is , is . As and , we find for the initial linear growth as long as the intermediate zone at is still present.
2.10. Initial Shape of Shock Front
We may now calculate the (stronger) shape of the stable linear displacing shock front that is expected to retain its shape as long as the intermediate layer and the supernatant layer are present, e.g. as long as shock’s rate of displacement is constant. In the shock, between and , the porosity is only a function of , thus in our pde is only a function of u, e.g. (u): Hence and . Thus, our pde is then only a function of u: . Integration once: , where c1 is an integration constant. But we know that for that or , so . Integrating once more for , as long as the intermediate zone and a supernatant zone are present, where and .
As the wave travels with constant shape, we may shift to pass at T = 0 the value = 1 at ξ = 1 and follow the curve: . The wave shape scales inversely proportional to δ. The flank is steeper when δ has a larger absolute value. For a sedimentation process at = 0.3 and δ = −30, we obtain Table 1 and Figure 2.
So, for , and when .
These parameter choices will be shown later to represent a virtual “sedimentation”, e.g. inverse of the actual experimentally observed creaming, of a Hellmann’s mayonnaise sample structured with whole egg (WE 6.2%) in a centrifuge at 1000
Table 1. The constant shape of the intermediate/supernatant boundary in linear displacement for = 0.3.
Figure 2. Shape of moving shock front at δ = −30 and = 0.3.
rpm. For this small sample rotating at large distance from the axis centripedal force can be considered constant, hence sedimentation and creaming are mirror images when density differences only differ in sign.
The slope of the wave front is indeed very steep over a large range of higher values, and there is a sharp nick towards the supernatant. Later, when the middle layer has disappeared, the porosity drops at the shock front and : hence or . In sedimentation or
The height of the clear supernatant layer Ω is easily monitored in the product and is often used to experimentally measure compaction over time in earth’s gravity or centrifuge. In this subsequent slowdown from elastic counter forces, we expect a displacement characterized by from the similarity with the diffusion and heat transport equations. We may anticipate this from a simplified model.
2.11. Simplified Model in Earth’s Gravity and/or at Small Density Differences in Limit δ ® 0
For small δ, like in earth’s gravity and or at small density differences, our pde is . For small changes of , we have approximately . The major effect in the small earth’s gravity is therefore that a value of will displace as square root of time , as observed in many other physical phenomena, like chemical diffusion, heat transport, et cetera.
For the slow 1-dimensional poroelastic drainage in the earth’s gravity field there is a functional dependence or scaling with a dimensionless group . This is an important scaling number for counter-current compaction drainage processes in the earth’s gravity in poroelastic media and is the equivalent of the Fourier number for heat transport.
A typical time for approach to equilibrium would be a displacement over height H, e.g. or a typical dimensionless time T = 0.43 when = 0.3. Suppose that the coefficient of consolidation is indeed a constant during the drainage process. Hence, the permeability reduction and the elasticity enhancement are assumed to be compensated by their multiplication. We take a fluid-filled poroelastic medium (like a gel, a porous or fibrous network, representing an emulsion and/or a particulate) that is compacted by drainage from one side (like pressing with an open filter that allows easy passage of fluid and does not allow any passage of the particulates). We may then solve the pde analytically. In the present case, we assume that: the initial condition is for ; boundary conditions for t>0 and for . It is therefore assumed that the material drains from one side in a one-dimensional way and that the material is so thick, or times are so short, that the drainage disturbance does not reach the back of the fabric, within the time span of observation. Standard mathematics leads to a solution . The amount of fluid flowing is . Liquid flows in the direction of lowest pressure and highest porosity. To make liquid flow out of the porous medium through z = 0, the pressure on the outside must be lower than inside. Then . Mass balance gives or . The liquid continues to flow with time at a diminishing rate. This is the case because it was assumed that the disturbance has not (yet) reached the other side of the medium. In our case, the production of fluid at z = 0 leads to a compaction shortening (e.g. depleted layer) of length , where Ω the compaction length. This simple model applies in processes like filtration compaction, machine and hand wash, et cetera. In this simplified gravity compaction model, the compaction shortening Ω is proportional to during the entire time span that the compaction disturbance has not reached the other side.
3.1. Height of Supernatant during Creaming of a Fibrous Nata de Coco Gel Network
The supernatant height Ω(t) = h(t) of an aqueous fibrous system of micro fibrillated bacterial cellulose (BC, Nata de Coco) and Carboxy Methyl Cellulose (CMC) at ratio 14/1, with captured soy bean oil (SB) droplets was followed in earth’s gravity as a function of time  , see Figure 3.
The volume of the sample is 20 ml, height ca. 42 mm. The experimental data of the initial rate as function of BC and Oil concentration are shown in Table 2;
The cake compacts proportionally to the square root of time:
The linear slopes are given in Table 3.
Figure 3. Creaming of Nata de Coco in water  .
Table 2. Initial rate as function of BC concentration  .
Figure 4. Supernatant height Ω of a creaming fibrous network of BC/CMC with SB oil as function of time  .
Figure 5. Supernatant height Ω of a creaming fibrous network of BC/CMC with SB oil as function of square root of time  .
Table 3. Slope as function of BC/CMC concentration  .
3.2. Creaming of Mayonnaise
Dressings with different formulations and varying processing were tested in a centrifuge to accelerate gravitational instabilities (sedimentation and creaming). These data are needed to predict the shelf life and to test structurants that allow the development of cheaper, lower carbon footprint, clean label, or better sustainable dressings. The LumiFuge centrifuge accelerates the study of gravitational instability behaviour of over a year into just 8 hours, due to a factor 1100 magnified gravitational force used in the centrifuge compared to earth’s gravity. Such data is often difficult to obtain in real time as the mayonnaises/dressings might be attacked from bacteria and/or fungi, might become chemically instable, or might change consistency on a smaller time scale, not directly related to gravitational instability during shelf life, but obscuring the results.
As the centrifuge data for creaming of the dressings is rather complex, including changing rate over time, changing rate with speed of revolution and with composition, it is difficult to pinpoint centrifuge conditions (time/speed) that can be directly related to the magnitude of the actual slow creaming in the earth’s gravitational field during shelf life. Therefore, a more elaborate analysis is warranted.
Several experiments focussed on the creaming of a Hellmann’s Real mayonnaise formulation WE 6.2%, both in earth’s gravity and in the LUMiFuge at 1000, 2000, 3000, and 4000 rpm; physical characterisation as follows.
3.2.1. Initial Porosity of WE 6.2%
The 67.5 wt% bean oil in the WE 6.2wt% formulation has a density of 0.9192 at 20˚C. Fresh whole egg consists typically of 74 wt% water, 11.8 wt% lipids and 12.8 wt% protein, and 1.4 wt% dissolved material (carbohydrate and minerals). A fresh egg has therefore 75.4 wt% water phase and 24.6 wt% solids. Some 34% of egg weight is the yolk. A fresh egg has a density of 1.033  . Therefore, 100 g of a WE 6.2wt% formulation consists of 6.2 g WE 6.2%, and contains 1.5 g eggs solids phase of density 1.132 g/cm3 and 4.7 g egg water of density 1.00. The rest of the WE 6.2% dressing is water plus dissolved solids = 100 − 67.5 − 6.2 = 26.3 wt%, with density 1.00. If we assume 100 g formulation, we have 67.5 g bean oil and 1.5 g as egg dispersed material and 26.3 + 4.7 = 31 g water phase. The volume of the dispersed solid is 67.5/0.9192 + 1.5/1.132 = 73.43 + 1.325 = 74.76 cm3 solids and 31/1.00cm3 water phase, thus in total 105.76 cm3. The initial volume percentage of dispersed material is therefore 74.76/105.76 = 70.7 vol%, hence 1 − = 0.707, and the volume percentage of water phase is 29.3, hence the original homogeneously mixed WE 6.2% formulation has an initial porosity of = 0.293.
3.2.2. Densities of WE 6.2%
The density of the aqueous fluid phase is around 1000 kg/m3, and the bean oil is 919.2 kg/m3. The density difference is therefore . We neglect the influence of the 1.5 wt% egg solids on the density difference, because the contribution is small and the egg white and yolk solids may redistribute between water and oil.
3.2.3. Viscosity of Fluid in WE 6.2%
The viscosity of the aqueous fluid phase is close to water, e.g. μ = 1 mPas.
3.2.4. Gravitation and LumiFuge
In our LUMiFuge centrifuge experiments, RGF can be 6-4000, vessel height H typically 22 mm, and axis of gyration typically 131 mm at the vessel’s base, or on average . We may assume that the LumiFuge’s centrifugal force is constant over the short height, which simplifies the mathematical analysis. When we assume an average radius of raver = 0.12 m and 3000 rpm leads to an average , where , then RGF = 1207 for WE 6.2% at 3000 rpm.
The LUMiFuge (Figure 6) provides the objective classification and quantification of de-mixing phenomena and fast determination of stability and shelf life of dispersions. Based on the patented STEP-Technology, the device determines simultaneously the de-mixing processes of 8 different dilute or concentrated suspensions/emulsions. The multi-sample analytical centrifuge LUMiFuge accelerates the de-mixing process up to 2300 times compared to traditional test tube tests. Stability tests and shelf life determinations of original dispersions are up to 2500 times faster than performed in a test tube under earth’s gravity.
The LUMiFuge (type 110, year of build 2004, LUM GmbH) is an analytical centrifuge with an opto-electronic sensor system that measures near infrared NIR transmission profiles along horizontally inserted sample tubes  . Data is acquired according to the patented STEP (Space and Time resolved Extinction Profiles) procedure  . Flat transparent sample tubes (Figure 7, LUM 2 mm PC Rectangular Synthetic Cell 110-131XX; 0.5 ml) are filled with 0.4 ml of sample, placed in the rotor and 225 NIR transmission profiles are recorded for each
Figure 6. The LUMiFuge centrifuge.
run. The tube’s base has a radius of about 131 mm, where the top of the evaluation range is near 104 mm radius. See for schematics Figure 8.
The LUMiFuge® software—SEPView V.6. (L.U.M. Ltd.)—calculates the space- and time-resolved transmission profiles. Standard front tracking analysis and instability index are tools to characterize sedimentation stability by the software, see for instance  . Creaming of a few prototype dressing formulations is followed in the LUMiFuge at 1000, 2000, 3000, and 4000 rpm (experiments done in duplicate or even triplicate) as well as by observation in earth’s gravity during shelf life over a month. An example of the LUMiFuge centrifuge data for Hellmann’s Real mayonnaise sample WE 6.2% is in Figure 9:
Each red curve is a transmission trace at a certain time, starting with the small nick curve on the right near 130 mm, progressing towards the green band near 125 mm in 2 hours. The steep slopes of the red (and ultimately green) lines are the top of the cleared (supernatant) zone building from the bottom of the tube in time. The creaming height Ω is (arbitrarily) taken at the distance between the curves at 10% transmission (indicated by the red arrow) and the bottom of the tube (starting with practically zero creaming height at time zero). The sample height is approximately 22 mm. The top is the stationary flank near 108 mm position, where all red and green curves superimpose, and is the position of the
Figure 7. The sample cells.
Figure 8. The LUMiFuge schematics. Our unit reaches 4000 rpm and 2300 ggrav at the base of the sample.
meniscus, at the top surface of the sample. A typical collection of results for samples running at 4000 rpm for 24 h is given in Figure 10.
Here we see the clear layer at the bottom of the sample. In some cases, here in the samples in the middle, we see a thin oil rim at the top of the sample. As these samples have been running for a long time (24 h), the clear layer has practically reached the maximum thickness. The clear basal layer indicates that even the tiniest particles are not only captured in the network but also move concurrently within it. When the water phase is not clear, the particles move more independently,
Figure 9. A sample of the LUMiFuge optical transmission versus position for Hellmann’s Real WE 6.2 wt% at 3000 rpm during 2 hours.
Figure 10. Samples that have run for 24 h at 4000 rpm. These samples showed the impact of varying structurants in some typical pourable dressings.
with small particles left behind: for example, the sample on the left. Such cases are not described by our modeling.
Creaming is also observed in the earth’s gravitational field with new samples of Hellmann’s Real in test tubes shown in Figure 11.
3.2.6. Experimental Data LumiFuge WE 6.2% at 3000 rpm
In these centrifuge experiments, the measured creaming height is sometimes indicated by h, but it actually always the height of the depleted layer Ω (in creaming) of our previous modeling. The creaming data of WE 6.2% at 3000 rpm from Figure 9 is shown in Figure 12.
For WE 6.2% at 3000 rpm, we see a linear slope ,
1 ml = 50/8 mm height. Initial height in gravity drainage of 10 ml (62.5 mm).
Figure 11. Test tubes used to measure creaming height in earth’s gravity.
Figure 12. The creaming height Ω as a radial position of the top of the water layer as function of time for WE 6.2% at 3000 rpm, .
e.g. 6.351 μm/s. From this initial slope we estimate k0 = 13.27 × 10−15 m2 (=13.27 mD), using Δρ = −80.8 kg/m3 and g = 11,843 m/s2, and μ = 1 mPas.
The graph shows the linear part (right-hand side) and the reduction of speed when the creaming zone reaches the top of the sample (left-hand side). The linear part allows us to calculate the creaming speed.
3.2.7. Experimental Data LumiFuge WE 6.2% at Several Speeds
The experimental LumiFuge data for WE 6.2% at 1000, 2000, 3000, and 4000 rpm are shown in Figure 15:
We see that all curves have essentially the same slope . There is a scatter in the waiting time. probably due to differences in spin up timing of centrifuge. The square root linear parts of are
Figure 13. The creaming height Ω of WE 6.2% as function of square root of time for 3000 rpm followed during 2 hours (experiments in duplicate).
Figure 14. Fit of linear part of creaming height Ω of WE 6.2% as function of square root of time for 3000 rpm (experiments in duplicate).
Our data are thus averaged:
3.2.8. Experimental Data WE 6.2% in Earth’s Gravity Field
We may compare these data with the visual observation of the creaming height
Figure 15. Log/log plot of WE 6.2% centrifuge data.
Table 4. WE 6.2% The fitted diffusive part of centrifuge data.
Table 5. WE 6.2% The fitted diffusive part of centrifuge data.
The first data points are tricky as the tube had a V-shaped base (first few mm), so the early rates are difficult to read.
3.2.9. Extrapolation Centrifuge Data of WE 6.2% to Reveal Shelf Life Stability
Combining the gravity and centrifuge rates versus square root of time in one graph, for WE 6.2%, we obtain Figure 17:
This graph shows that we cannot use a linear extrapolation to predict the creaming in the earth’s gravity field. But if we plot the data on a log/log plot, see Figure 18, we find surprisingly:
So, for WE 6.2% . In SI units
with r2 = 0.978 and sdev = 0.084. At each experiment, there is a definite scaling of Ω2/t proportional to RGF which becomes identical with earth’s gravity data at RGF = 1. When RGF increases, we would expect the experiment to go faster as we exert a larger force ω2R. It is therefore not unlikely that Ω2/t is some function of RGF, e.g. even appearing to be directly proportional to RGF. We expect similar responses for other structurants than whole egg at 6.2 wt%, with a different proportionality constant.
This equation already allows us to predict the (square root of time) rate in earth’s gravity (RGF = 1) from the rates measured in centrifuge by a simple
Table 6. Creaming of WE 6.2% in earth’s gravity field.
Figure 16. Creaming position of WE 6.2% in earth’s gravity field.
Figure 17. Creaming slope position of WE 6.2% as function of RGF, the g-force relative to gravity acceleration (e.g. 9.812 m/s2): .
Figure 18. Log/log plot creaming slope of WE 6.2% as function of .
log/log extrapolation, e.g. from a linear fit on a log/log plot of Ω2/t versus RGF.
For WE 6.2% at 3000 rpm, the linear part in square root of time quantified to and 247.12 , on average 258 in these two repeat experiments. We use the average slope at 3000 rpm of , e.g. . The diffusive part is approximately , hence for the diffusive/dispersive front . Using the known values , k0 = 13.273 mD and μ = 1 mPas, we estimate ε0 = 3870 N/m2. Thus, for creaming of WE 6.2% at 3000 rpm, δ = 261.
We may look at the equilibrium Θ of WE 6.2% at 3000 rpm: , where is the ultimate porosity at the entrance (ξ = 0). So, for δ = 261 (β = 5.44) and Θequil ~ 5/22 = 0.2273 we find ~ 0.106. This is not far from the observed value.
Using a stepwise integration of with small steps Δξ = 0.01 from ξ = 0 and = = 0.107, we may calculate the equilibrium curves. This is done for δ = −261, δ = −30 (and = 0.3) in Figure 19 and Figure 20:
We may use k0 = 13.7 × 10−15 m2 and ε0 = 3870 N/m2 based on data at 3000 rpm to predict the δ parameters for other gravities, as is done in Table 7:
Figure 19. The predicted equilibrium curve for WE 6.2% at 3000 rpm.
Likewise we find for 1000 rpm:
Figure 20. The predicted equilibrium curve for WE 6.2% at 1000 rpm.
Table 7. Creaming Hellmann’s Real mayonnaise WE 6.2% at equilibrium prediction based on experimental data at 3000 rpm.
NB. In sedimentation, Δρ is positive and thus δ is just the negative value, and hence the porosity profile is the vertical mirror image, when the acceleration (buoyancy) can be considered constant (earth’s gravity field) or approximated to be constant (small sample in centrifuge). We may conclude that for the mayonnaise centrifuge experiments with small H and Δρ and thus, small value of δ, the diffuse square root dependence dominates quickly in the response.
3.2.10. Waiting Time
In all experiments, we see some time has passed before the cleared layer becomes visible, e.g. the waiting time. Part of that is the short time needed to rise the porosity to unity in the top of the vessel. In our experiments another part is the spin up time of the centrifuge. If a structured sample is slowly shaken just before the experiment to even out the porosity (or set up vertically after first slowly rolling horizontally), it might take some time before the firmed structure is developed, which might contribute to the waiting time.
The waiting time can be prolonged in the experiment by adding hygroscopic ingredients, like starch, that adsorb the initial free water, by swelling the porous matrix. This adsorption technique is often applied to extend shelf life in mayonnaises. The swelling suppresses the build-up of a depleted layer until the diminishing rate of swelling is equal to the depletion rate, from which moment on, the depleted layer starts to build, quickly accelerating to a constant rate, as swelling is decaying exponentially in time. The swelling and prolonged waiting time is not in the model yet.
3.3. Subsidence in the Groningen Gasfield, The Netherlands
When we look at, for instance, subsidence in the earth’s gravity field, like Groningen’s gasfield  , as shown in Figure 21:
We see that the subsidence response Ω/t even after several years of constant gas production was still in its initial linear constant slope transient. Nowadays, the gas production has declined, so a new, smaller linear trend slope of constant Ω/t must have started, that will transform eventually in a dispersive trend Ω2/t before the equilibrium subsidence depth is approached after many years. An independent determined value of the average elasticity of the rock formation
Figure 21. Subsidence data for Groningen’s gasfield.
might be used to predict this later time Ω2/t approach to equilibrium in the future and thus predict the decay of subsidence in the far future.
4. Numerical Simulation
In our numerical simulations, we describe the centrifuge WE 6.2% creaming experiments in the equivalent “sedimentation” mode by simulating with the opposite (negative) value of δ, which means for the small samples in the centrifuge a simple reversal of the direction of the vertical axis, e.g. simulating the mirror image in vertical direction.
Porosity condition , spatial vertical height coordinate , time 0 ≤ T. All symbols dimensionless. Initial condition (t = 0) for WE6% mayonnaise constant, e.g. = 0.293, boundary conditions no flow over boundary and , or when < 1. As a required Lipschitz continuous initial non-decreasing profile in sedimentation, the sedimentation process starts with a constant porosity, as conforms to the experiments. The solution is first constant and later in sedimentation either constant or increasing, producing an upward curve from compacting cake on bottom (left) side to the top (right) side of the graph below. The porosity always increases quickly towards = 1 at the top of the vessel, with a steep wave slope at the boundary. The solution, when porosity has reached unity in the top of the vessel, continues moving in a sharp boundary front (Stefan condition) towards supernatant layer with = 1. The front moves inward while the supernatant layer grows in height. A weak solution given as Rankine-Hugoniot conditions over the sharp front where and the values at the top of the cake. Earlier, the front is towards the intermediate layer, hence is , e.g. a front displacing with a constant rate and supernatant thickness growing linearly. When the height of the intermediate layer has become zero, the supernatant grows gradually slower in height until the equilibrium porosity profile is reached, with a balance between gravity and elastic forces. Typical values for δ in “sedimentation” of WE 6% mayonnaise are the negative of those given in Table 7 for several speeds. Scaling height H = 0.022 m, z = Hξ; scaling time τ = 45,477 s, t = Tτ.
The first rise of porosity in the end plane from = to = 1 over a short time can be generated as a strong solution (no jumps yet) in space fixed coordinates. The time involved is natural waiting time. An example calculated with COMSOL5.3a for δ = −30, = 0.3, approximately equivalent to WE 6.2% at 1000 rpm, is shown in Figure 22. Four time-steps of 0.006 to dimensionless time T = 0.024 have just passed the = 1 level in the exit plane, indicating that for this system the waiting time is just above T = 0.018:
After the waiting time, the further build-up of the cake in the presence of an intermediate constant layer can be generated in space fixed coordinates as a strong solution in a virtually infinite vessel by assuming at the exit . This works until the loss of porosity in the cake compaction is equal to the gain of porosity in the supernatant, e.g. the point in time where the intermediate layer is eroded. During this erosion period, the growth of the supernatant thickness Θ at porosity over length is linear in time, as generated by the Rankine-Hugoniot jump condition , e.g. . An example calculated with COMSOL5.3a for δ = −30, = 0.3, approximately equivalent to WE 6.2%
Figure 22. Comsol simulation of waiting time for δ = −30, = 0.3 (equivalent to centrifuge “sedimentation” experiment WE 6.2% 1000 rpm). No flow through boundaries. Four time-steps shown.
at 1000 rpm is shown in Figure 23. Ten time-steps of 0.01 to dimensionless time T = 0.1 are generated in a long vessel ξ = 5:
We see that at T = 0.01 (end of first time-step), the cake has already reached the normal exit boundary at ξ = 1, indicating that in a normal vessel length ξ = 1 the intermediate layer has eroded at T = 0.01 for WE 6.2% at 1000 rpm and the cake had already reached the thin supernatant layer. Hence, T = 0.01 is near the end of the linear growth period of the supernatant.
The subsequent rate of cake formation slows, initially proportional to square root of time, and must probably be generated in material coordinates as a weak solution with a Stefan boundary condition. We tried a strong solution, but that is not stable in the supernatant layer. By mere luck we generated a solution that was reasonably stable for longer time. This example, calculated with COMSOL5.3a for δ = −30, = 0.3, approximately equivalent to WE 6.2% at 1000 rpm is shown in Figure 24. Ten time-steps of 0.005 to dimensionless time T = 0.05 remain reasonably stable:
Here we see a build-up of the supernatant layer and a diminishing of the supernatant depth growth rate at later time (T > 0.01). The supernatant porosities are fluctuating ( ¹ 1), and not reaching = 1, with a steep slope and a sharp nick to steady = 1 in the supernatant, indicating numerical instability. The later time displacements, and other simulation attempts, were not accurate due to the numerical errors in the supernatant range.
During slow down, improving the accuracy of the simulations might confirm the expected Θ2 µ T, e.g. Ω2 µ t behaviour. Because the displacement profile up to equilibrium has some shape similarity, we might try a self-similar transformation to predict larger time behaviour. This is not yet done.
Figure 23. Comsol simulation interim δ = −30, = 0.3 (equivalent to centrifuge experiment WE 6.2% 1000 rpm) Left bc no flow, right bc is Dirichlet = 0.3. Ten time-steps shown.
Figure 24. Comsol simulation δ = -30, = 0.3 (equivalent to centrifuge experiment WE 6.2% 1000 rpm) No flow through boundaries. Ten time-steps shown.
Transformation to Material Coordinates
A trick might be used to avoid the difficult Stefan moving boundary problem by stretching the caked layer to remove the supernatant layer from the domain of integration. We have a dispersed material volume balance such that any two material-bound spatial coordinates at a small distance dξ0 compact as
. We might therefore define new spatial coordinates
that would be invariant under the compaction: our pde is then for
, 0 < T:
, with boundary conditions (bc) at χ = 0
= 1 at χ = 1, and initial condition T = 0,
for 0 < χ < 1. Expressed in coordinate χ, the zone with dispersed material is virtually stretched, such that the compaction is compensated to keep the total length of that zone constant. By rescaling
, such that
, 0 < T, with
such that bc u = 0 is
and bc at
= 1. Over domain (0
Notably, at equilibrium , e.g. integrated . At the entrance = for u = 0 or which can be solved for . For = 0.3 and δ = −30 we find = 0.179, which is indeed close to the predicted equilibrium curve (based on an estimated value of Θequil). We may now calculate Θequil accurately by using and to find , so , thus, at equilibrium: . For = 0.179, = 0.3, δ = −30, we find ξflank = 0.70, e.g. Θequil = 0.30. Hence, we expect the experimentally estimated (long time) “equilibrium” thickness of about 5.5 mm to grow to 6.6 mm at full equilibrium (H was 22 mm). The earlier estimate Θequil = 0.25 was probably too short. There is no need to estimate that equilibrium supernatant length from the experiment, because it is determined by the 4 independent Pi parameters and can thus, together with the equilibrium intercept , be directly calculated from those parameters.
Using solution profile (u, T) at any given T, we may re-tabulate the proper curve (ξ, T) by the transformation . It is important to realize that we were indeed able to calculate all characteristics of the compaction curves from waiting, transient time, up to equilibrium, with only the 4 independent dimensionless groups (t/τ, z/H, δ, and ) as given by the Pi theorem. The parameters themselves were calculated using only one value of the seven physical parameters [ , Δρ, g, H, k, ε, μ] that appear in determining those 4 dimensionless groups. The only two difficult physical parameters to determine independently are permeability k and porosity ε. However, these can be determined using the experimental data, from the early linear slope of Θ vs t and later linear slope of Θ vs respectively.
The simulation data in material height coordinates u for sedimentation WE 6.2% at several rpm’s are given in Figure 25. These simulations are reasonably accurate but are also time consuming (several hours on a 64 bit W10 AMD FX-8150 Eight-Core 8 Gb ram 3.6 GHz bench top computer). An overview of several representative simulations (COMSOL5.3a):
Apparently, the flank and slope become steeper and flatter for increasing , with both the value and trend in the porosity of the intercept at u = 0 agreeing with the experimental observations. As the vertical differences in subsequent time-steps rapidly diminish, at an order of half of the difference in the previous time-step, the equilibrium curve is expected to be lower, approximately less than or equal to the shift during the last recorded time-step, hence already very close to the last and lowest recorded curve. As noticed before, the proper curves (ξ, T) can be found by the transformation at any T. An example of this conversion is given in Figure 26 for WE6.2% at 3000 rpm, δ = −260 at
Figure 25. Simulations δ = −10, −30, −60, and −260 at = 0.3 of cake porosity profiles in u coordinates as function of dimensionless time in steps ΔT. Here horizontal coordinate is material height coordinate u.
Figure 26. The porosity (ξ,T) for δ = −260 and = 0.3 for 10 steps ΔT = 0.005 to T = 0.05.
= 0.3 for T from 0 to 0.05 in steps of ΔT = 0.005:
This gives an indirect way to resolve (ξ,T) without stepping into the intricate Stefan moving boundary problem, and still find a rather strong solution. The physics of the moving boundary problem keeps the jump sharp even in these rather strong solutions, as was indeed observed in many experiments. NB. The initial flat profile was modified by a small upward curve near the vessel’s top, to suppress oscillations and blow-up at the start.
5. Bead and Spring Model
The poroelastic counter-current compaction can be seen as a porous elastic weak “Bead and Spring” model, subjected to a gravity field, with the coordinate system (e.g. coordinate z) attached to the dispersed medium (the bead and spring system) of lateral area A, as shown in Figure 27.
Rotate graph 90 degrees counter-clockwise for actual position in sedimentation, 90 degrees clockwise for creaming. The porous medium in the vessel is assumed to consist of N cells filled with a (repeating bead and spring) elastic porous medium. Initially, all cell “boundaries” are equidistant at distance Δz0 = H/N and porosity . The left side is the bottom wall attached to the first cell. The right-hand side N is initially a bead next to the other vessel wall, without being attached to the wall (or simply ending in a water-filled medium to the right of that point N). The clearing between the right-hand side wall and the right-hand side bead is the thickness Ω of the depletion layer (initially zero).
When gravity is applied, the spring chain will fall (displace to the left on the graph) compressing the elastic springs (near the bottom wall) and pushing the fluid
Figure 27. Bead and Spring Model in natural dimensions and scaled dimensionless.
inside the dispersed system to the right, displacing the coordinate zN, when the dispersed material of the last cell moves to the left (starting with the same fall displacement rate as all other beads). The volume in each cell is initially where with A area, H length of vessel, and N number of beads and springs, with dispersed volume . The fluid movement will displace cell j boundaries zj and zj−1, such that the same volume of dispersed material stays in between the boundaries in the cell j, that has at any time a volume of . The coordinate zj is attached to the dispersed material and moves co-current with that dispersed material at the upper boundary of cell j. The porosity of cell j is thus at any time
. The volume of dispersed material in a cell is conserved, hence the volume change in a cell is the volume of fluid that is pushed out of the cell. The difference of these flows through the cell boundaries compresses the cell fluid volume (per unit time), hence per unit area: . Thus and of course the solid is moving counter-currently: . Darcy has the factor 1/2 in this counter-current flow: . Using Carman-Kozeny en Van Wyk and the dimensionless scaling as used before: . Therefore, our scaled pde has not changed: . Rearranging is . We may therefore use a volume-conservative numerical scheme: for a time-step ΔT between Tk and Tk+1 for each spring element j = 1, N where ξ0 = 0 and initially springs equidistant zj = j/N. This procedure is often used in literature for this Richards’ type of equations  .
The bead and spring model is a weak solution with an infinitely sharp flank towards the supernatant. Therefore, the Rankine-Hugenoit jump condition at the flank towards the supernatant (with conditions and = 1) applies: . Here, and are taken at the top of the cake, e.g. the low porosity side of the vertical flank (the upstream side) and ξ* a coordinate attached to the material in the dispersed phase.
We will try a different route to solve the pde using the weak bead and spring model. We had derived above that for a coordinate zj attached to the dispersed material moving co-current with dispersed material at the upper cell wall of cell j that or expressed in our scaled coordinates . The fluid flow through a horizontal level in the cake must displace the dispersed material over a distance −dξ*. For a spring element therefore, , where and and ΔT has a time-step k to k + 1. The flow rate is calculated by . The average porosity in the spring element j is given by . If assuming that the flow rate changes approximately linear over each cell, the average porosity represents the middle value of the cell, e.g. . We need to interpolate these middle cell porosities at to the value at the upstream side cell boundary j, e.g. at the position of “the bead” at . For cell wall j: . For the lower first cell boundary j = 0: . Extrapolating for the upper cell upper wall N: . The last grid point ξ*(N) coincides continually with the start of the depletion layer at ξ = 1 − Θ in the original space fixed grid. The volume conservation of the fluid can be checked by , which always equals , whereby is non-negative and Θ increases in each step, otherwise, the simulation stops. An improvement might be to impose the required monotony in and , and require the , the maximum slope at equilibrium. This is not implemented yet. Even when starting with equidistant spatial mesh, the springs near the bottom wall compact quickly to create a non-equidistant mesh. The stable size of the time-step can be determined by the displacement rate . For the cake with , a maximum rate at or occurs. We expect sufficient small displacements to maintain sufficient stability if ΔT = 0.5Δξ/Umax. A Fortran simulation for the bead and spring model for sedimentation WE 6.2% in centrifuge for = 0.293 for δ = −0.22, −29, −115, −261 and −464, for N = 10 or 20, appears reproducibly
Figure 28. Dimensionless scaled supernatant height Q(δ,N) = W/H where N is the number of dimensionless spatial steps (N = 10, 20) in the simulation as function of dimensionless scaled time T.
consistent up to a reasonable time, as shown in Figure 28:
The experimental data are scaled dimensionless by τ = 45,477 s and H = 0.022 m. In comparing the simulation with the experimental data at 3000 rpm, it appears that most misfit is due to the uncertainty of the zero time: a small waiting time is expected. The observed order of magnitude of 100 seconds seems rather large but not unrealistic considering the preparation and design of the centrifuge experiment where the centrifuge has to gain speed. The fit can be improved by adapting the roughly guessed simulation parameters τ, δ, and incorporating a waiting time, but this was not pursued further, since the reproducibility of the experiments is low (see the experimental section, Figure 15). The detailed simulated profile for 3000 rpm is shown for δ = −261, N = 10 in Figure 29:
The elementary bead and spring model is fast, but accuracy is uncertain, especially at longer time.
The 1-dimensional consolidation theory in soil is often based on Terzaghi’s theory that assumes that the soil is fully saturated with water, that water and soil are incompressible, the effective stress is the sum of the total stress in the rock and the pore pressure, strains in soil are small, Darcy’s Law applies to hydraulic flow, permeability and volume compressibility are constant, and there is a unique relation between porosity and stress, independent of time. Our theory can cope with nonlinear changing of permeability and elasticity during the consolidation processes over larger, non-linear strains and stresses, and correctly incorporates the buoyancy force.
Figure 29. The porosity development of WE 6.2% in centrifuge (δ = −261 = 0.2983).
Even the Terzaghi theory is better than expected from the small amplitude linearization. When written as , the coefficient of consolidation is a product of kε which is a better constant over a wider range of consolidation due to the compensating effects in the porosity dependence of the product of k and ε.
The same pde equation applies in washing, when a falling cloth filled with washing liquid hits the rotating drum, creating on impact a pressure 1/2 ρv2 from its fall speed v, or an article is squeezed, or hit by a stick, exerting again a pressure, such that the fluid in the pores moves suddenly with respect to the dispersed material that then compacts and builds an elastic counter force that slows the drainage. The initial fast drainage relative to the fibres gives the best wash action, dislodging dirt from the fibres. In partially aerated clothes, e.g. with foaming, the wash action on impact repeats—less efficiently—in subsequent patches of fluid. This explains why fully saturated clothes wash less efficiently. Not only is foam in washing liquor aesthetically pleasing, but the foam-inducing surfactant facilitates lowering of the interfacial tension which dislodges dirt, keeping it in suspension, and assists in the mechanics of the washing process by creating repeated air patches within the clothes.
Another action is the permeable piston applying pressure on a fluid-infiltrated dispersed system (cake), such that it compacts in a filtration operation, draining fluid from dispersed material. Such processes are described by the same equation, with the appropriate adapted initial and boundary conditions.
Even the reverse process, swelling of dehydrated cake, is a process that might be described by the same type of equations and relations, at least when the process is slow and reversible (elastic).
7. Impact and Outlook
There is a wide range of measures to stabilize dispersions/emulsions based on experience in the field that can now be understood and quantified according to our model: make fluid viscous and turn fluid into a gel during storage; add structurants to liquids (starch, alginate, micro-fibrillated cellulose, gums like LBG and guar); add a hygroscopic material (starch) to absorb the first free liquid; increase viscosity, decrease elasticity, permeability and density difference. We may now measure stability quickly in fast centrifuge and translate results straightforward into slow changes in shelf life in earth’s gravity.
There are a large collection of models and literature on sedimentation, creaming and synaeresis in porous, granular, or fibrous media, for granulates, for suspensions, for dispersions, or for emulsions. The solutions presented there often apply to a smaller range and are usually based on semi-empirical interaction parameters for particle-particle interactions between specific particles. Our model uses a trick to circumvent the specific disperse interaction and focus on the fluid flow in poroelastic media, using smart scaling relations that make the theory more universal and applicable in many different formulations. Phenomenological models for thickening during sedimentation of flocculated suspensions are generally based on the Kynch model, where the compaction reaches a fixed lower porosity limit in the cake to avoid further compaction     . In our model, that condition is reached by the local dynamic balance of buoyancy and elastic forces in the compaction drainage of porous media, using widely applicable smooth permeability and elasticity functions that at the same time determine the range of applicability, together with the simplifying assumptions mentioned in the introduction (§1.1). The counter-current flow without segregation requires initial porosities below say typically 45%, or small values of parameter δ. Discrete droplet aeration asks for an overrun below 50% - 150%, depending on level of air stabilizers. In soils, the residual gas saturations should be below typically 35% to trap the gas as discrete bubbles.
Extending the model by additional physics, like starting with a particular monotone profile, adding compressibility of fluid or solids, anisotropy, solving the full centrifuge equations in cylindrical symmetry, seems straightforward. Even replacing the Carman-Kozeny or Van Wijk equations respectively by the permeability and elasticity characterisations using other dedicated equations as a function of porosity can easily be done. The reverse process of drainage compaction, e.g. imbibition of fluid infiltrating into swelling porous media, is also important (in reservoir and civil engineering, like capturing rainfall), but is not modeled here, as plastic or elastic responses require different solutions. The reversible elastic mode of dynamic counter-current imbibition will resemble the reverse of our counter-current drainage compaction solutions. The extension to more than one flowing (connected) fluid phase (aqueous, and/or oleic, and/or gaseous) requires extra physics (capillarity, hysteresis, relative permeability) and additional models (Van Genuchten, Corey, Brooks). In our model, additional fluids are only present as disconnected dispersed bubbles (overrun) or droplets (emulsions) that are trapped in the dispersed phase and flow co-current with that dispersed phase, and counter-current with the draining fluid.
I would like to thank Z. Salazar, E. M. Hutter and K. P. Velikov for providing the experimental data and for the stimulating discussion.
 Barry, S.I. and Mercer, G.B. (1999) Flow and Deformation in Poroelasticity: I Unusual Exact Solutions. Mathematical and Computer Modelling, 30, 23-29.
 Henderson, N., Brêttas, J.C. and Sacco, W.F. (2010) A Three-Parameter Kozeny–Carman Generalized Equation for Fractal Porous Media. Chemical Engineering Science, 65, 4432–4442.
 Kuentz, M. and Roethlisberger, D. (2003) Rapid Assessment of Sedimentation Stability in Dispersions Using Near Infrared Transmission Measurements during Centrifugation and Oscillatory Rheology. European Journal of Pharmaceutics and Biopharmaceutics, 56, 355-361.
 Zeng, X. and Decker, M. (2009) Improving the Numerical Solution of Soil Moisture-Based Richards Equation for Land Models with a Deep or Shallow Water Table. Journal of Hydrometeorology, 10, 308-319.
 Bürger, R. and Concha, F. (1998) Mathematical Model and Numerical Simulation of the Settling of Flocculated Suspensions. International Journal of Multiphase Flow, 24, 1005-1023.