ABSTRACT

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 compaction in counter-current fluid drainage. Processes captured are amongst others sedimentation, creaming and subsidence. The most important limiting prerequisite is that the incompressible dispersed medium is sufficiently structured and/or concentrated that it compacts during slow drainage, without segregation in sizes or in components. 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.

Keywords:

Poroelastic, Compaction, Consolidation, Gravity, Centrifuge, Sedimentation, Creaming, Subsidence, Emulsion, Dispersion

1. Introduction

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.

1.1. Assumptions

The continuous incompressible fluid (aqueous or oleic) is assumed Newtonian.

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

1.2. Method

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.

1.3. Nomenclature

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 [1] [2] [3] [4] [5] . 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 $\frac{\partial \varphi }{\partial t}+\frac{\partial q}{\partial z}=0$ 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 $q=-\frac{k}{2\mu }\frac{\text{d}P}{\text{d}z}$ with $\frac{\text{d}P}{\text{d}z}=\Delta \rho \left(1-\varphi \right)g+\epsilon \frac{\partial \varphi }{\partial z}$, where k is permeability of the porous medium, μ the viscosity of the fluid, Δρ the averaged density difference between the dispersed solids and the fluid, $\varphi =\varphi \left(t,z\right)$ 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 $\frac{\partial \varphi }{\partial t}=-\frac{\partial }{\partial z}\left[-\frac{k}{2\mu }\left\{\Delta \rho \left(1-\varphi \right)g+\epsilon \frac{\partial \varphi }{\partial z}\right\}\right]$.

2.1. Darcy Equation

The Darcy equation was derived originally for fluid flow in stagnant soil or rock as $q=-\frac{k}{\mu }\frac{\text{d}P}{\text{d}z}$. 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 $2q=-\frac{k}{\mu }\frac{\text{d}P}{\text{d}z}$. This explains the factor 1/2 in front of Darcy in counter-current flow; e.g. $q=-\frac{k}{2\mu }\frac{\text{d}P}{\text{d}z}$. 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: $\Delta \rho =\frac{{\displaystyle {\sum }_i\Delta {\rho }_{i,s}{\varphi }_{i,s}}}{{\displaystyle {\sum }_i{\varphi }_{i,s}}}$, which is a constant during the compaction. Here, $\Delta {\rho }_{i,s}$ 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 ${\varphi }_{s,i}$ 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 $k=\frac{2{\varphi }^3}{{\alpha }_c{a}^2{\left(1-\varphi \right)}^2}$, where a is the grain surface area per unit of grain volume and a_c an experimental correlation coefficient [6] [7] . 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 $k=\frac{c{\varphi }^3}{a_g^2}\approx k_0{\left(\frac{\varphi }{{\varphi }_0}\right)}^3{\left(\frac{1-{\varphi }_0}{1-\varphi }\right)}^2$, where k₀ is the permeability at a reference porosity ${\varphi }_0$. The elasticity of a fibrous network pack varies with compression. To describe the elasticity, we will use the Van Wyk’s heuristic equation [8] . For a fibrous fabric, it relates the applied pressure P to the inverse cube of the volume V, the intrinsic uncompressed volume V₀₀, corrected for the incompressible (high pressure solid limit) volume V_s as follows [9] [10] : $P=K_{c}E\frac{m^3}{{\rho }^3}\left(\frac{1}{{\left(V-V_s\right)}^3}-\frac{1}{{\left(V_{00}-V_s\right)}^3}\right)$. 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 K_c characterising the fibres, of order 0.01 that will vary with fibre orientation. Now, the overall elasticity of the gel network is defined by $\epsilon =-V\left(\frac{\partial P}{\partial V}\right)$, thus using the modified Van Wyk’s equation $\epsilon =3V{K}_{c}E\frac{m^3}{{\rho }^3}\left(\frac{1}{{\left(V-V_s\right)}_{}^4}\right)$. The fluid-filled pore space is $V-V_s=V\varphi$. In vertical compression, the volume of particulates is conserved; $V_s=V\left(1-\varphi \right)=V_{00}\left(1-{\varphi }_{00}\right)$. Thus, $\epsilon =3K_{c}E\frac{m^3}{{\rho }^3}\left(\frac{1}{V^3}\frac{1}{{\varphi }^4}\right)=3K_{c}E\frac{m^3}{V_{00}^3{\rho }^3}{\left(\frac{1-\varphi }{1-{\varphi }_{00}}\right)}^3\frac{1}{{\varphi }^4}$. So, between two states in compaction $\epsilon ={\epsilon }_0{\left(\frac{1-\varphi }{1-{\varphi }_0}\right)}^3{\left(\frac{{\varphi }_0}{\varphi }\right)}^4$. We selected the reference state as the porosity condition, where we can measure the value of porosity dependent physical parameters like k₀ and ε₀. 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 $\frac{\partial \varphi }{\partial t}=\frac{k_0{\epsilon }_0}{2\mu }\frac{\partial }{\partial z}\left[\frac{k}{k_0}\left\{\frac{\Delta \rho g}{{\epsilon }_0}\left(1-\varphi \right)+\left(\frac{\epsilon }{{\epsilon }_0}\right)\frac{\partial \varphi }{\partial z}\right\}\right]$, or using Carman-Kozeny and Van Wyk: $\frac{\partial \varphi }{\partial t}=\frac{k_0{\epsilon }_0}{2\mu }\frac{1}{H^2}\frac{\partial }{\partial \frac{z}{H}}\left[{\left(\frac{\varphi }{{\varphi }_0}\right)}^3{\left(\frac{1-{\varphi }_0}{1-\varphi }\right)}^2\left\{\frac{\Delta \rho gH}{{\epsilon }_0}\left(1-\varphi \right)+{\left(\frac{1-\varphi }{1-{\varphi }_0}\right)}^3{\left(\frac{{\varphi }_0}{\varphi }\right)}^4\frac{\partial \varphi }{\partial \frac{z}{H}}\right\}\right]$, where we introduced the height of the aquifer or vessel H as spatial scaling. If we now introduce a dimensionless scaling parameter $\delta =\frac{{\left(1-{\varphi }_0\right)}^3}{{\varphi }_0^4}\frac{\Delta \rho gH}{{\epsilon }_0}$, e.g. negative in sedimentation, a dimensionless time $T=\frac{t}{\tau }=\frac{k_0{\epsilon }_0}{2\mu }\frac{1}{H^2}\frac{{\varphi }_0}{1-{\varphi }_0}t$, and a dimensionless scaled vertical spatial coordinate $\xi =\frac{z}{H}$ with $0\le \xi \le 1$, then pde becomes fully dimensionless as $\frac{\partial \varphi }{\partial T}=\frac{\partial }{\partial \xi }\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}$. If we define a dimensionless fluid flow rate $Q^{\prime }$ as $-Q^{\prime }=\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}$, the dimensionless volume balance becomes $\frac{\partial \varphi }{\partial T}=-\frac{\partial Q^{\prime }}{\partial \xi }$, where $Q^{\prime }=q\frac{\tau }{H}$. 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 $\varphi$, we have 7 independent parameters [ ${\varphi }_0$, Δρ, 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 ${\Pi }_1={\varphi }_0$,${\Pi }_2=\xi =z/H$,${\Pi }_3=T=t/\tau$ and ${\Pi }_4=\delta$ determine the dependent variable $\varphi$. Because ${\varphi }_0$ is dimensionless, it is a Pi number in itself. As z and H have the same dimension, their ratio is a Pi number [11] . Thus, we need only one consistent set of the 7 independent parameter values [ ${\varphi }_0$, Δρ, 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: $F^2=\frac{\rho L^2{v}^2}{\rho L^{3}g}=\frac{v^2}{Lg}$. In the same way, the Cauchy number is the ratio of the inertia force over the elastic force $C^2=\frac{\rho L^2{v}^2}{E{L}^2}=\frac{\rho v^2}{E}$. 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: $\frac{C^2}{F^2}=\frac{\text{gravitational}\text{force}}{\text{elastic}\text{force}}=\frac{\rho L^{3}g}{E{L}^2}=\frac{\rho gL}{E}$, or in our case more specifically $\frac{\Delta \rho gH}{{\epsilon }_0}$. This is proportional to our $\delta =\frac{\Delta \rho gH}{{\epsilon }_0}\frac{{\left(1-{\varphi }_0\right)}^3}{{\varphi }_0^4}$. The δ is therefore also a measure of the ratio of gravitational and elastic force. For ${\varphi }_0$ = 0.3 we expect a balance when δ = 4.2.

2.5. Initial and Boundary Conditions

The initial condition is (for now) constant porosity ${\varphi }_0$ 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: ${\displaystyle \underset{0}{\overset{1}{\int }}\varphi \text{d}\xi }={\varphi }_0$. 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. $\frac{\partial \varphi }{\partial T}\le 0$, which means that in view of the volume conservation $\frac{\partial \varphi }{\partial T}=-\frac{\partial Q^{\prime }}{\partial \xi }$ that $\frac{\partial Q^{\prime }}{\partial \xi }\ge 0$ at that cake boundary. If the pde prevails up to the boundary, and $\varphi$ can be differentiated twice, then at that no-flow boundary ${Q^{\prime }}_{bc}=0={\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}}_{bc}$, or ${\left\{\frac{\partial \varphi }{\partial \xi }\right\}}_{bc}={\left\{-\frac{\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}\right\}}_{bc}$. 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 ${\varphi }_0$ towards $\varphi$ ® 1, so there is $\frac{\partial \varphi }{\partial T}\ge 0$ and thus $\frac{\partial Q^{\prime }}{\partial \xi }\le 0$. Thus, as long as the porosity function is not saturated ( $\varphi$ < 1), the no-flow boundary condition at the depletion side is ${Q^{\prime }}_{bc}=0={\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}}_{bc}$ and thus ${\left\{\frac{\partial \varphi }{\partial \xi }\right\}}_{bc}={\left\{\frac{-\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}\right\}}_{bc}$. 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 ( $0<\varphi <1$ ), the porosity obeys the pde over the entire range, hence ${\displaystyle \underset{0}{\overset{1}{\int }}\frac{\partial \varphi }{\partial T}\text{d}\xi }=0={\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}}_{\xi =0}^{\xi =1}=-\left[Q^{\prime }\left(\xi =1\right)-Q^{\prime }\left(\xi =0\right)\right]=0$, e.g. $Q^{\prime }\left(\xi =0\right)=Q^{\prime }\left(\xi =1\right)=0$, confirming the Robin nonlinear boundary conditions at both sides.

However, as the porosity reaches $\varphi$ = 1 at an infinite thin depletion/supernatant layer, the boundary condition is a no-flow condition $Q^{\prime }=0$ at the bottom flank of the depletion layer of thickness Θ, where $\varphi$ = 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 $\varphi$ = 1 constantly. The lower boundary of the depleted layer is an almost vertical moving front ending in a sharp nick towards $\varphi$ = 1. We may interpret this moving boundary as a space shifting of a Dirichlet boundary condition $\varphi$ = 1 (over a spatial range with an extend Θ from the top of the vessel).

If we write the pde as $\frac{\partial \varphi }{\partial T}+\frac{\partial }{\partial \xi }\left\{\varphi U\left(\varphi \right)\right\}=\frac{\partial }{\partial \xi }\left\{D\left(\varphi \right)\frac{\partial \varphi }{\partial \xi }\right\}$, 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 $U\left(\varphi \right)=\frac{-\delta {\varphi }^2}{1-\varphi }$ is an increasing function of $\varphi$. The LHS leads therefore to a shock front (kinematic wave) propagating with a speed U( ${\varphi }_0$ ); the sharp front separates the homogeneous suspension of initial porosity ${\varphi }_0$ 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 $D\left(\varphi \right)=\left(\frac{1-\varphi }{\varphi }\right)$ becomes infinitely small for high values of $\varphi$ 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 ${\left\{\frac{\partial \varphi }{\partial \xi }\right\}}_{bc}={\left\{\frac{-\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}\right\}}_{bc}$ that the flank of the cleared zone is also very steep for high values of $\varphi$. 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 [12] . 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 ${\varphi }_0$ and for a given ratio δ of buoyancy and elastic force. Repeating the calculations for other values of ${\varphi }_0$ 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 $\varphi$ (T, ξ) determined by a combination ( ${\varphi }_0$, δ), there is in real space and time a fan of curves $\varphi$ (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 $\varphi$ = ${\varphi }_0$ 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 $q=-\frac{k}{2\mu }\frac{\text{d}P}{\text{d}r}$ with $\frac{\text{d}P}{\text{d}r}=\Delta \rho \left(1-\varphi \right)g+\epsilon \frac{\partial \varphi }{\partial r}$. In the gravitational field, the acceleration $g=g{\hat{e}}_z$, is practically constant (g = −9.8125 m/s²), but in a centrifuge the centripetal acceleration is a function of r: $g=g{\hat{e}}_r$, with $g={\omega }^{2}r$, 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: $\frac{\partial \varphi }{\partial t}=-\frac{1}{r}\frac{\partial \left(rq\right)}{\partial r}$. For a small rectangular vessel at larger distance from the axis (say 4 times its length), the volume balance equations can be approximated as $\frac{\partial \varphi }{\partial t}=-\frac{1}{r}\frac{\partial \left(rq\right)}{\partial r}\approx -\frac{1}{\bar{r}}\frac{\partial \left(\bar{r}q\right)}{\partial r}=-\frac{\partial q}{\partial r}=\frac{\partial q}{\partial Z}$, where $\bar{r}$ 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 $Z=r-r_0$. In sedimentation, r₀ is the outer radius of the sample and $r_0=\bar{r}+H/2$. In creaming, r₀ is the inner radius of the sample and $r_0=\bar{r}-H/2$ and $r=r_0+Z$, with Δρ being negative. So, for such a small vessel spinning at large distances from the rotation axis $\frac{\partial \varphi }{\partial t}=\frac{\partial q_w}{\partial Z}$ and $\frac{\text{d}P}{\text{d}r}=-\frac{\text{d}P}{\text{d}Z}=\Delta \rho \left(1-\varphi \right){\omega }^2\left(r_0-Z\right)-\epsilon \frac{\partial \varphi }{\partial Z}\approx \Delta \rho \left(1-\varphi \right){\omega }^2\bar{r}-\epsilon \frac{\partial \varphi }{\partial Z}$. So $-{\omega }^2\bar{r}$ has for the small sample in the centrifuge the same role as g_grav, but often represents a much larger gravity in fast spinning centrifuges: ${\omega }^2\bar{r}=-RGF\cdot g_{grav}$ where RGF a large positive number, or in such centrifuges: $\frac{\text{d}P}{\text{d}Z}\approx \Delta \rho \left(1-\varphi \right)RGF\cdot g_{grav}+\epsilon \frac{\partial \varphi }{\partial Z}$. 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, τ, ${\varphi }_0$, and δ.

2.7. Equilibrium Profile

In the ultimate equilibrium, the fluid flow $Q^{\prime }=0$ over the entire vessel, hence $\frac{\text{d}\varphi }{\text{d}\xi }=\frac{-\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}$. Integration yields the equilibrium profile $\delta \left(\xi -{\xi }_k\right)=-\left[-\frac{1}{3{\varphi }^3}+\frac{1}{{\varphi }^2}-\frac{1}{\varphi }+\frac{1}{3{\varphi }_k^3}-\frac{1}{{\varphi }_k^2}+\frac{1}{{\varphi }_k}\right]$, where ${\varphi }_k$ some reference point. If we take the reference point at the point $\xi =1-{\Theta }_{equil}$ where $\varphi \to \text{1}$ in equilibrium: $\xi -\left(1-{\Theta }_{equil}\right)=-\frac{1}{\delta }\left[-\frac{1}{3{\varphi }^3}+\frac{1}{{\varphi }^2}-\frac{1}{\varphi }+\frac{1}{3}\right]$, so at the bottom ξ = 0 where $\varphi ={\varphi }_{0\infty }$ : $1-{\Theta }_{equil}=-\frac{1}{\delta }\left[-\frac{1}{3}+\frac{1}{3{\varphi }_{0\infty }^3}-\frac{1}{{\varphi }_{0\infty }^2}+\frac{1}{{\varphi }_{0\infty }}\right]$. The supernatant is a clear fluid layer with $\varphi$ = 1 over scaled distance ${\Theta }_{equil}$.

In radial symmetry in equilibrium

$\frac{\partial \varphi }{\partial t}=-\frac{1}{r}\frac{\partial \left(rq\right)}{\partial r}=-\frac{1}{r}\frac{\partial }{\partial r}\left(-r\frac{k}{2\mu }\frac{\text{d}P}{\text{d}r}\right)=\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{k}{2\mu }\left(\Delta \rho \left(1-\varphi \right)g(r)+\epsilon \frac{\partial \varphi }{\partial r}\right)\right)=0$ leads to $\Delta \rho \left(1-\varphi \right)g\left(r\right)+\epsilon \frac{\text{d}\varphi }{\text{d}r}=\text{const}$. At the axis r = 0, the porosity is (virtually) maintained at the initial value $\varphi ={\varphi }_0$ as there is still the initial pressure irrespective of rotation speed, and ${\left(\frac{\text{d}\varphi }{\text{d}r}\right)}_{r=0}\equiv 0$ from the axial symmetry, hence const = 0, which leads to the same condition as in constant gravity, but now z is replaced by r: $\Delta \rho \left(1-\varphi \right)g(r)=-\epsilon \frac{\text{d}\varphi }{\text{d}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: $\frac{\frac{1}{2}{r}^2}{H\bar{r}}=-\frac{1}{\delta }\left[-\frac{1}{3{\varphi }^3}+\frac{1}{{\varphi }^2}-\frac{1}{\varphi }+\frac{1}{3{\varphi }_0^3}-\frac{1}{{\varphi }_0^2}+\frac{1}{{\varphi }_0}\right]$. 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 $\varphi$ to be continuous and at least twice differentiable over its domain. The solution of our pde: $\frac{\partial \varphi }{\partial T}=\frac{\partial }{\partial \xi }\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\partial \varphi }{\partial \xi }\right)=-\frac{\partial Q^{\prime }}{\partial \xi }$ 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: $\frac{\partial }{\partial T}{\displaystyle \underset{{\xi }_1}{\overset{{\xi }_2}{\int }}\varphi \text{d}\xi }={\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\partial \varphi }{\partial \xi }\right)}_{{\xi }_1}^{{\xi }_2}=-\left(Q^{\prime }\left({\xi }_2,T\right)-Q^{\prime }\left({\xi }_1,T\right)\right)$. This integral describes the same physics but needs to be only once differentiable. This is a weak formulation and is approximate. $\frac{\partial }{\partial T}{\displaystyle \underset{0}{\overset{\xi }{\int }}\varphi \text{d}\xi }={\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\partial \varphi }{\partial \xi }\right)}_0^{\xi }=-Q^{\prime }\left(\xi ,T\right)$. 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 $\frac{\partial \varphi }{\partial T}=-\frac{\partial Q^{\prime }}{\partial \xi }$ can be rewritten as a weak shock front speed ${\left(\frac{\text{d}\xi }{\text{d}T}\right)}_{\varphi }=-\frac{\Delta Q^{\prime }}{\Delta \varphi }$, where $-Q^{\prime }=\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}$. Over a shock front the volume balance leads to Rankine-Hugoniot condition $\frac{\text{d}\xi }{\text{d}T}=-\frac{{Q^{\prime }}_{\downarrow }-{Q^{\prime }}_{\uparrow }}{{\varphi }_{\downarrow }-{\varphi }_{\uparrow }}$. For the initial linear displacement cleared front Θ between ${\varphi }_0$ and depleted layer with 1 with respectively $Q^{\prime }={Q^{\prime }}_0=\frac{-\delta {\varphi }_0^3}{\left(1-{\varphi }_0\right)}$ and $Q^{\prime }=0$, we find $\frac{\text{d}\xi }{\text{d}T}=-\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}$. Since the dimensionless clear zone length in sedimentation is $\Theta =1-\xi$, is $\frac{\text{d}\Theta }{\text{d}T}=\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}=\frac{\Delta \rho gH}{{\epsilon }_0}\frac{\left(1-{\varphi }_0\right)}{{\varphi }_0}$. As $T=\frac{t}{\tau }=\frac{k_0{\epsilon }_0}{2\mu }\frac{1}{H^2}\frac{{\varphi }_0}{1-{\varphi }_0}t$ and $\Omega =H\Theta$, we find for the initial linear growth $\frac{\text{d}\Omega }{\text{d}t}=\frac{1}{2}\Delta \rho g\frac{k_0}{\mu }$ as long as the intermediate zone at $\varphi ={\varphi }_0$ 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 $U=\frac{\text{d}\xi }{\text{d}T}=\frac{-\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}$ is constant. In the shock, between $\varphi ={\varphi }_0$ and $\varphi =\text{1}$, the porosity $\varphi \left(\xi \right)=\varphi \left(\xi +UT\right)$ is only a function of $u=\xi +UT$, thus $\varphi$ in our pde $\frac{\partial \varphi }{\partial T}=\frac{\partial }{\partial \xi }\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\partial \varphi }{\partial \xi }\right)=-\frac{\partial Q^{\prime }}{\partial \xi }$ is only a function of u, e.g. $\varphi$ (u): Hence $\frac{\partial \varphi }{\partial T}=\frac{\partial \varphi }{\partial u}\frac{\partial u}{\partial T}=U\frac{\text{d}\varphi }{\text{d}u}$ and $\frac{\partial \varphi }{\partial \xi }=\frac{\partial \varphi }{\partial u}\frac{\partial u}{\partial \xi }=\frac{\text{d}\varphi }{\text{d}u}$. Thus, our pde is then only a function of u: $U\frac{\text{d}\varphi }{\text{d}u}=\frac{\text{d}}{\text{d}u}\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\text{d}\varphi }{\text{d}u}\right)$. Integration once: $U\varphi =\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\text{d}\varphi }{\text{d}u}\right)+c_1=-Q^{\prime }+c_1$, where c₁ is an integration constant. But we know that for $\varphi ={\varphi }_0$ that $-Q^{\prime }=\frac{-\delta {\varphi }_0^3}{1-{\varphi }_0}$ or $-\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}{\varphi }_0=\frac{\delta {\varphi }_0^3}{1-{\varphi }_0}+c_1$, so $-\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}\varphi =\left(\frac{\delta {\varphi }^3}{1-\varphi }+\frac{1-\varphi }{\varphi }\frac{\text{d}\varphi }{\text{d}u}\right)-\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}$. Integrating once more $u=\xi +UT={\displaystyle \int \frac{\frac{1-\varphi }{\varphi }}{-\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}\varphi -\frac{\delta {\varphi }^3}{1-\varphi }+\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}}\text{d}\varphi }$ for ${\varphi }_0<\varphi <\text{1}$, as long as the intermediate zone and a supernatant zone are present, where $u=\xi +\frac{-\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}T$ and $\begin{array}{l}u=\frac{{\left(1-{\varphi }_0\right)}^2}{\delta {\varphi }_0^3}\{-\frac{2{\varphi }_0\text{arctanh}\left[\frac{2\varphi {\left(1-{\varphi }_0\right)}^2+{\varphi }_0-2{\varphi }_0^2}{{\varphi }_0\sqrt{4{\varphi }_0-3}}\right]}{\left({\varphi }_0-3\right)\sqrt{4{\varphi }_0-3}}+\ln \varphi \\ \begin{array}{c}\\ \\ \\ \end{array}-\frac{\left({\varphi }_0-1\right)\ln \left(\varphi -{\varphi }_0\right)}{{\varphi }_0-3}+\frac{\ln \left[{\varphi }^2{\left(1-{\varphi }_0\right)}^2+{\varphi }_0^2+\varphi \left({\varphi }_0-2{\varphi }_0^2\right)\right]}{{\varphi }_0-3}\}+c_2\end{array}$.

As the wave travels with constant shape, we may shift to pass at T = 0 the value $\varphi$ = 1 at ξ = 1 and follow the curve: $\left(1-\xi \left(\varphi \right)\right)=\frac{1}{\delta }{\displaystyle \underset{\varphi }{\overset{1}{\int }}\frac{\frac{1-\varphi }{\varphi }}{\frac{{\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}\left(1-\varphi \right)-\frac{{\varphi }^3}{1-\varphi }}\text{d}\varphi }={\displaystyle \underset{\varphi }{\overset{1}{\int }}\frac{1}{\frac{{\varphi }_0^3\varphi }{{\left(1-{\varphi }_0\right)}^2}-\frac{{\varphi }^4}{{\left(1-\varphi \right)}^2}}\text{d}\varphi }$. The wave shape scales inversely proportional to δ. The flank is steeper when δ has a larger absolute value. For a sedimentation process at ${\varphi }_0$ = 0.3 and δ = −30, we obtain Table 1 and Figure 2.

So, for $\delta =-\text{3}0$,$\xi \left(\varphi =0.50\right)=0.988$ and $\xi \left(\varphi =0.301\right)=0.34$ when $\xi \left(\varphi =1\right)=1$.

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

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 $\varphi$ values, and there is a sharp nick towards the supernatant. Later, when the middle layer has disappeared, the porosity drops at the shock front $\varphi <{\varphi }_0$ and $\frac{\partial \varphi }{\partial \xi }>0$ : hence ${\left(\frac{\text{d}\xi }{\text{d}T}\right)}_{\varphi }=-\frac{\Delta Q^{\prime }}{\Delta \varphi }$ or ${\left(\frac{\text{d}\xi }{\text{d}T}\right)}_{\varphi }=-\frac{{Q^{\prime }}_{\downarrow }-{Q^{\prime }}_{\uparrow }}{{\varphi }_{\downarrow }-{\varphi }_{\uparrow }}=-\frac{0-\left(\frac{-\delta {\varphi }^3}{1-\varphi }-\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right)}{1-\varphi }=\frac{-\delta {\varphi }^3}{{\left(1-\varphi \right)}^2}-\left(\frac{1}{\varphi }\right)\frac{\partial \varphi }{\partial \xi }$. In sedimentation ${\left(\frac{\text{d}\Theta }{\text{d}T}\right)}_{\varphi }=-\frac{\text{d}\xi }{\text{d}T}=\frac{\delta {\varphi }^3}{{\left(1-\varphi \right)}^2}+\left(\frac{1}{\varphi }\right)\frac{\partial \varphi }{\partial \xi }$ or $\frac{1}{H}{\left(\frac{\text{d}\Omega }{\text{d}t}\right)}_{\varphi }=\frac{k_0{\epsilon }_0}{2\mu }\frac{1}{H^2}\frac{{\varphi }_0}{1-{\varphi }_0}{\left[\frac{{\left(1-{\varphi }_0\right)}^3}{{\varphi }_0^4}\frac{\Delta \rho gH}{{\epsilon }_0}\frac{{\varphi }^3}{{\left(1-\varphi \right)}^2}+\left(\frac{1}{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right]}_{\varphi =\varphi \left(\text{at}\text{low}\text{side}\text{of}\text{shock}\right)}$

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 $\Omega ~\sqrt{t}$ 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 $\frac{\partial \varphi }{\partial t}\approx \frac{k_0{\epsilon }_0}{2\mu }\frac{\partial }{\partial z}\left\{\left(\frac{1-\varphi }{1-{\varphi }_0}\right)\left(\frac{{\varphi }_0}{\varphi }\right)\frac{\partial \varphi }{\partial z}\right\}$. For small changes of $\varphi$, we have approximately $\frac{\partial \varphi }{\partial t}\approx \frac{k_0{\epsilon }_0}{2\mu }\frac{{\partial }^2\varphi }{\partial z^2}$. The major effect in the small earth’s gravity is therefore that a value of $\varphi$ will displace as square root of time $z~\sqrt{t}$, 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 $\frac{k_0{\epsilon }_0}{2\mu }\frac{t}{z^2}$. 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. $t\sim {\frac{2\mu H}{k_0{\epsilon }_0}}^2$ or a typical dimensionless time T = 0.43 when ${\varphi }_0$ = 0.3. Suppose that the coefficient of consolidation $D_s=\frac{k_0{\epsilon }_0}{2\mu }$ 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 $\varphi \left(z,0\right)={\varphi }_0$ for $0<z<\infty$ ; boundary conditions $\varphi \left(0,t\right)=0$ for t>0 and $\partial \varphi /\partial z\to 0$ for $z\to \infty$. 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 $\varphi \left(z,t\right)={\varphi }_{0}erf\left(\frac{z}{2\sqrt{D_{s}t}}\right)$. The amount of fluid flowing is $q\left(z,t\right)=-\frac{k_0}{2\mu }\frac{\partial P}{\partial z}=-D_s\frac{\partial \varphi }{\partial z}$. 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 $q\left(z,t\right)=-D_s\frac{\partial \varphi }{\partial z}=-D_s\frac{\partial }{\partial z}\left({\varphi }_{0}erf\left(\frac{z}{2\sqrt{D_{s}t}}\right)\right)$. Mass balance gives ${\displaystyle \int q}\left(z=0,t\right)\text{d}t=-{\displaystyle \underset{0}{\overset{t}{\int }}{D}_s{\varphi }_0{\left[\frac{2{\text{e}}^{-\frac{z^2}{4D_{s}t}}}{\sqrt{\text{π}}}\frac{1}{2\sqrt{D_{s}t}}\right]}_{z=0}\text{d}t}$ or $\frac{Q\left(z=0,t\right)}{A}={\displaystyle \int q\left(z=0\right)\text{d}t}=-\frac{2}{\sqrt{\text{π}}}\sqrt{D_s}{\varphi }_0\sqrt{t}$. 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 ${\displaystyle \int q}\left(z=0,t\right)\text{d}t=\Omega \left(t\right)=-\frac{2}{\sqrt{\text{π}}}\sqrt{D_s}{\varphi }_0\sqrt{t}$, 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 $\sqrt{t}$ during the entire time span that the compaction disturbance has not reached the other side.

3. Experiments

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 [13] , 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 creaming data for 0.06 wt% BC/CMC with 1 wt% oil is reproduced in Figure 4 and Figure 5 (experiments in triplicate):

The cake compacts proportionally to the square root of time:

The linear slopes are given in Table 3.

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 [14] . 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/cm³ 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 cm³ solids and 31/1.00cm³ water phase, thus in total 105.76 cm³. The initial volume percentage of dispersed material is therefore 74.76/105.76 = 70.7 vol%, hence 1 − ${\varphi }_0$ = 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 ${\varphi }_0$ = 0.293.

3.2.2. Densities of WE 6.2%

The density of the aqueous fluid phase is around 1000 kg/m³, and the bean oil is 919.2 kg/m³. The density difference is therefore $\Delta \rho =\left(1-0.9192\right)1000=80.8\text{kg}/{\text{m}}^{\text{3}}$. 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 $\bar{r}=120\text{mm}$. We may assume that the LumiFuge’s centrifugal force $g={\omega }^2\bar{r}$ is constant over the short height, which simplifies the mathematical analysis. When we assume an average radius of r_aver = 0.12 m and 3000 rpm leads to an average $g=11843=1207\times \left(-g_{grav}\right)\text{m}/{\text{s}}^{\text{2}}$ _, where $g_{grav}=-9.8125\text{m}/{\text{s}}^{\text{2}}$, then RGF = 1207 for WE 6.2% at 3000 rpm.

3.2.5. LumiFuge

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 [15] . Data is acquired according to the patented STEP (Space and Time resolved Extinction Profiles) procedure [16] . 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

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 [16] . 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

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,

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 $-\frac{\text{d}R}{\text{d}t}=\frac{\text{d}\Omega }{\text{d}t}=16462\left(\text{mm}/\text{month}\right)$,

e.g. 6.351 μm/s. From this initial slope $\frac{\text{d}\Omega }{\text{d}t}=-\frac{1}{2}\frac{k_0\Delta \rho g}{\mu }$ we estimate k₀ = 13.27 × 10⁻¹⁵ m² (=13.27 mD), using Δρ = −80.8 kg/m³ and g = 11,843 m/s², and μ = 1 mPas.

The two traces of the creaming height as function of the square root of time for WE 6.2% at 3000 rpm are given in Figure 13 and Figure 14 (duplicate experiments, axes reversed).

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 ${\xi }^2~\Delta \rho g_{grav}\text{RGF}T$. There is a scatter in the waiting time. probably due to differences in spin up timing of centrifuge. The square root linear parts of $\xi ~\sqrt{\Delta \rho g_{grav}\text{RGF}T}$ are

fitted, where $\text{RGF}=\left|\frac{{\omega }^{2}R}{g_{grav}}\right|$. Data is collected in Table 4 and averaged in Table 5.

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

Ω in a test tube as function of storage time (in earth’s gravitational field), see Table 6 and Figure 16.

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% $\frac{{\Omega }^2}{t}\approx \left(48.3{\text{mm}}^{\text{2}}/\text{month}\right)\cdot \text{RGF}$. In SI units

${}^{10}\mathrm{l}og\left[\Omega /\sqrt{t}/\left(\text{m}/{\text{s}}^{\text{0}\text{.5}}\right)\right]=\left(-5.36\pm 0.08\right)+\left(0.528\pm 0.03\right){}^{10}\mathrm{l}og\left(\text{RGF}\right)$ with r² = 0.978 and sdev = 0.084. At each experiment, there is a definite scaling of Ω²/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 ω²R. It is therefore not unlikely that Ω²/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

log/log extrapolation, e.g. from a linear fit on a log/log plot of Ω²/t versus RGF.

For WE 6.2% at 3000 rpm, the linear part in square root of time quantified to $\text{d}h/\text{d}\sqrt{t}=268.96$ and 247.12 $\text{mm}/\sqrt{\text{month}}$, on average 258 $\text{mm}/\sqrt{\text{month}}$ in these two repeat experiments. We use the average slope at 3000 rpm of $\Omega /\sqrt{t}=258\text{mm}/\sqrt{\text{month}}$, e.g. ${\Omega }^2/t=2.568\times {10}^{-8}{\text{m}}^2/\text{s}$. The diffusive part is approximately $\frac{\partial \varphi }{\partial t}\approx \frac{k_0{\epsilon }_0}{2\mu }\frac{{\partial }^2\varphi }{\partial z^2}$, hence for the diffusive/dispersive front $\frac{z^2}{t}=\frac{{\Omega }^2}{t}\approx \frac{k_0{\epsilon }_0}{2\mu }$. Using the known values ${\Omega }^2/t=2.57\times {10}^{-8}{\text{m}}^2/\text{s}$, k₀ = 13.273 mD and μ = 1 mPas, we estimate ε₀ = 3870 N/m². Thus, for creaming of WE 6.2% at 3000 rpm, δ = 261.

We may look at the equilibrium Θ of WE 6.2% at 3000 rpm: $1-{\Theta }_{equil}=-\frac{1}{\delta }\left[-\frac{1}{3}+\frac{1}{3{\varphi }_{0\infty }^3}-\frac{1}{{\varphi }_{0\infty }^2}+\frac{1}{{\varphi }_{0\infty }}\right]$, where ${\varphi }_{0\infty }$ is the ultimate porosity at the entrance (ξ = 0). So, for δ = 261 (β = 5.44) and Θ_equil ~ 5/22 = 0.2273 we find ${\varphi }_{0\infty }$ ~ 0.106. This is not far from the observed value.

Using a stepwise integration of $\Delta \varphi =\frac{-\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}\Delta \xi$ with small steps Δξ = 0.01 from ξ = 0 and $\varphi$ = ${\varphi }_{0\infty }$ = 0.107, we may calculate the equilibrium curves. This is done for δ = −261, δ = −30 (and ${\varphi }_0$ = 0.3) in Figure 19 and Figure 20:

We may use k₀ = 13.7 × 10⁻¹⁵ m² and ε₀ = 3870 N/m² based on data at 3000 rpm to predict the δ parameters for other gravities, as is done in Table 7:

Likewise we find for 1000 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 [17] , 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 Ω²/t before the equilibrium subsidence depth is approached after many years. An independent determined value of the average elasticity of the rock formation

might be used to predict this later time Ω²/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 $0\le \varphi \le \text{1}$, spatial vertical height coordinate $0\le \xi \le \text{1}$, time 0 ≤ T. All symbols dimensionless. Initial condition (t = 0) for WE6% mayonnaise $\varphi$ constant, e.g. ${\varphi }_0$ = 0.293, boundary conditions no flow over boundary $Q^{\prime }\left(\xi =0\right)=0$ and $Q^{\prime }\left(\xi =1\right)=0$, or ${\left(\frac{\partial \varphi }{\partial \xi }\right)}_{bc}=\frac{-\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}$ when $\varphi$ < 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 $\varphi$ = 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 $\varphi$ = 1. The front moves inward while the supernatant layer grows in height. A weak solution given as Rankine-Hugoniot conditions over the sharp front ${\left(\frac{\text{d}\xi }{\text{d}T}\right)}_{\xi =1-\Theta }=\frac{Q^{\prime }}{1-\varphi }$ where $Q^{\prime }$ and $\varphi$ the values at the top of the cake. Earlier, the front is towards the intermediate layer, hence is ${\left(\frac{\text{d}\xi }{\text{d}T}\right)}_{\xi =1-\Theta }=\frac{{Q^{\prime }}_0}{1-{\varphi }_0}$, e.g. a front displacing with a constant rate and supernatant thickness $\left|\Theta \right|$ growing linearly. When the height of the intermediate layer has become zero, the supernatant grows gradually slower in height $\left|\Theta \right|$ 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 $\varphi$ = ${\varphi }_0$ to $\varphi$ = 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, ${\varphi }_0$ = 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 $\varphi$ = 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 $\varphi ={\varphi }_0$. 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 $\varphi =1$ over length $\left|\Theta \right|$ is linear in time, as generated by the Rankine-Hugoniot jump condition $\frac{\text{d}\Theta }{\text{d}T}=\frac{-{Q^{\prime }}_0}{1-{\varphi }_0}=\frac{\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}$, e.g. $\left|\Theta \right|=\frac{-\delta {\varphi }_0^3}{{\left(1-{\varphi }_0\right)}^2}T$. An example calculated with COMSOL5.3a for δ = −30, ${\varphi }_0$ = 0.3, approximately equivalent to WE 6.2%

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 $\left|\Theta \right|$ 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, ${\varphi }_0$ = 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 ( $\varphi$ ¹ 1), and not reaching $\varphi$ = 1, with a steep slope and a sharp nick to steady $\varphi$ = 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 Θ² µ T, e.g. Ω² µ 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.

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ξ₀ compact as $\text{d}\xi =\frac{1-{\varphi }_0}{1-\varphi }\text{d}{\xi }_0$. We might therefore define new spatial coordinates $\text{d}\chi =\frac{1-\varphi }{1-{\varphi }_0}\text{d}\xi$ that would be invariant under the compaction: our pde is then for $\varphi =\varphi \left(T,\chi \right)$, 0 < T: $\frac{\partial \varphi }{\partial T}=\frac{1-{\varphi }_0}{1-\varphi }\frac{\partial }{\partial \chi }\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-{\varphi }_0}{\varphi }\right)\frac{\partial \varphi }{\partial \chi }\right\}$, with boundary conditions (bc) at χ = 0 $\left(1-{\varphi }_0\right)\frac{\partial \varphi }{\partial \chi }=\frac{-\delta {\varphi }^4}{1-\varphi }$ and $\varphi$ = 1 at χ = 1, and initial condition T = 0, $\varphi$ = ${\varphi }_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 $u=\chi /\left(1-{\varphi }_0\right)$, such that $\varphi =\varphi \left(T,u\right)$, 0 < T, with $\frac{\partial \varphi }{\partial T}=\frac{1}{1-\varphi }\frac{\partial }{\partial u}\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1}{\varphi }\right)\frac{\partial \varphi }{\partial u}\right\}$ such that bc u = 0 is $\frac{\partial \varphi }{\partial u}=\frac{-\delta {\varphi }^4}{1-\varphi }$ and bc at $u=1/\left(1-{\varphi }_0\right)$ is $\varphi$ = 1. Over domain (0 0 < u < 1 / ( 1 − ϕ 0 ) , and Lipschitz condition $0\le \varphi \le \text{1}$.<∞),>

Notably, at equilibrium $\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1}{\varphi }\right)\frac{\partial \varphi }{\partial u}=0$, e.g. integrated $-\delta \left(\frac{1}{1-{\varphi }_0}-u\right)=\frac{1}{6}+\frac{1}{3{\varphi }^3}-\frac{1}{2{\varphi }^2}$. At the entrance $\varphi$ = ${\varphi }_{0\infty }^{}$ for u = 0 or $0<-\delta \left(\frac{1}{1-{\varphi }_0}\right)=\frac{1}{6}+\frac{1}{3{\varphi }_{0\infty }^3}-\frac{1}{2{\varphi }_{0\infty }^2}$ which can be solved for ${\varphi }_{0\infty }^{}$. For ${\varphi }_0$ = 0.3 and δ = −30 we find ${\varphi }_{0\infty }^{}$ = 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 $\text{d}u=\frac{\text{d}\chi }{1-{\varphi }_0}$ and $\text{d}\chi =\frac{1-\varphi }{1-{\varphi }_0}\text{d}\xi$ to find $\text{d}u=\frac{1-\varphi }{{\left(1-{\varphi }_0\right)}^2}\text{d}\xi$, so $\text{d}\xi =\frac{{\left(1-{\varphi }_0\right)}^2}{1-\varphi }\text{d}u=-\frac{{\left(1-{\varphi }_0\right)}^2}{1-\varphi }\frac{1-\varphi }{\delta {\varphi }^4}\text{d}\varphi =\frac{{\left(1-{\varphi }_0\right)}^2}{-\delta }\frac{1-\varphi }{{\varphi }^4}\text{d}\varphi$, thus, at equilibrium: ${\xi }_{flank}\left(\varphi =1\right)=1-{\Theta }_{equil}=\frac{{\left(1-{\varphi }_0\right)}^2}{-\delta }\left[\frac{1}{6}+\frac{1}{3{\varphi }_{0\infty }^3}-\frac{1}{2{\varphi }_{0\infty }^2}\right]$. For ${\varphi }_{0\infty }^{}$ = 0.179, ${\varphi }_0$ = 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 ${\varphi }_{0\infty }^{}$, be directly calculated from those parameters.

Using solution profile $\varphi$ (u, T) at any given T, we may re-tabulate the proper curve $\varphi$ (ξ, T) by the transformation $\text{d}\xi =\frac{{\left(1-{\varphi }_0\right)}^2}{1-\varphi }\text{d}u$. 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 ${\varphi }_0$ ) as given by the Pi theorem. The parameters themselves were calculated using only one value of the seven physical parameters [ ${\varphi }_0$, Δρ, 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 $\sqrt{t}$ 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 $\left|\delta \right|$, 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 $\varphi$ (ξ, T) can be found by the transformation $\text{d}\xi =\frac{{\left(1-{\varphi }_0\right)}^2}{1-\varphi }\text{d}u$ at any T. An example of this conversion is given in Figure 26 for WE6.2% at 3000 rpm, δ = −260 at

${\varphi }_0$ = 0.3 for T from 0 to 0.05 in steps of ΔT = 0.005:

This gives an indirect way to resolve $\varphi$ (ξ,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 Δz₀ = H/N and porosity ${\varphi }_0$. 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

inside the dispersed system to the right, displacing the coordinate z_N, 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 $V_{cell}=A\Delta z_0$ where $A\Delta z_0=AH/N$ with A area, H length of vessel, and N number of beads and springs, with dispersed volume $V_{cell,dispersed}=\left[1-{\varphi }_0\right]V_{cell}=\left[1-{\varphi }_0\right]A\Delta z_0$. The fluid movement will displace cell j boundaries z_j and z_j−1, such that the same volume of dispersed material $V_{cell,dispersed}=\left[1-{\varphi }_0\right]A\Delta z_0$ stays in between the boundaries in the cell j, that has at any time a volume of $V_{cell,j}=A\left(z_j-z_{j-1}\right)$. The coordinate z_j 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

${\varphi }_j=\frac{V_{cell,j}-V_{cell,dispersed}}{V_{cell,j}}=1-\frac{\left(1-{\varphi }_0\right)\Delta z_0}{z_j-z_{j-1}}$. 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: $q_{w,j}-q_{w,j-1}=-\frac{\text{d}{z}_j}{\text{d}t}+\frac{\text{d}{z}_{j-1}}{\text{d}t}$. Thus $q_{w,j}=-\frac{\text{d}{z}_j}{\text{d}t}$ and of course the solid is moving counter-currently: $q_{s,j}=+\frac{\text{d}{z}_j}{\text{d}t}$. Darcy has the factor 1/2 in this counter-current flow: $q_w=-\frac{\text{d}z}{\text{d}t}=-\frac{k}{2\mu }\frac{\text{d}{P}_w}{\text{d}z}=-\frac{1}{2}\frac{k}{\mu }\left(g\Delta \rho \left(1-\varphi \right)+\epsilon \frac{\text{d}\varphi }{\text{d}z}\right)$. Using Carman-Kozeny en Van Wyk and the dimensionless scaling as used before: $-Q^{\prime }=\frac{\text{d}\xi }{\text{d}T}=\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\text{d}\varphi }{\text{d}\xi }$. Therefore, our scaled pde has not changed: $\frac{\partial \varphi }{\partial T}=\frac{\partial }{\partial \xi }\left\{\frac{\delta {\varphi }^3}{1-\varphi }+\left(\frac{1-\varphi }{\varphi }\right)\frac{\partial \varphi }{\partial \xi }\right\}=-\frac{\partial Q^{\prime }}{\partial \xi }$. Rearranging is $\frac{\partial \xi }{\partial T}=-\frac{\partial Q^{\prime }}{\partial \varphi }$. We may therefore use a volume-conservative numerical scheme: $\left({\xi }_j-{\xi }_{j-1}\right)\frac{\left({\varphi }_{j-\frac{1}{2}}^{k+1}-{\varphi }_{j-\frac{1}{2}}^k\right)}{\Delta T}=-{\left(Q^{\prime }\right)}_j^{k+1}+{\left(Q^{\prime }\right)}_{j-1}^{k+1}$ for a time-step ΔT between T_k and T_k+1 for each spring element j = 1, N where ξ₀ = 0 and initially springs equidistant z_j = j/N. This procedure is often used in literature for this Richards’ type of equations [18] .

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 $Q^{\prime }=0$ and $\varphi$ = 1) applies: $\frac{\text{d}\Theta }{\text{d}T}=-{\left(\frac{\text{d}\xi }{\text{d}T}\right)}_{flank}=-{\left(\frac{\text{d}{\xi }^{*}}{\text{d}T}\right)}_{flank}=\frac{{Q^{\prime }}_{\downarrow }-{Q^{\prime }}_{\uparrow }}{{\varphi }_{\downarrow }-{\varphi }_{\uparrow }}={\left(\frac{-Q^{\prime }}{1-\varphi }\right)}_{upstream}$. Here, $Q^{\prime }$ and $\varphi$ 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 z_j attached to the dispersed material moving co-current with dispersed material at the upper cell wall of cell j that $q_{w,j}=-\frac{\text{d}{z}_j}{\text{d}t}$ or expressed in our scaled coordinates ${Q^{\prime }}_j=-\frac{\text{d}{\xi }_j^{*}}{\text{d}T}$. The fluid flow $Q^{\prime }\text{d}T$ through a horizontal level in the cake must displace the dispersed material over a distance −dξ^*. For a spring element therefore, $\Delta {\xi }_j^{*k+1}={\xi }_j^{*k+1}-{\xi }_{j-1}^{*k+1}=\Delta {\xi }_j^{*k}-\left({Q^{\prime }}_j-{Q^{\prime }}_{j-1}\right)\Delta T$, where ${\xi }_{j=0}^{*}=0$ and ${Q^{\prime }}_{j=0}=0$ and ΔT has a time-step k to k + 1. The flow rate ${Q^{\prime }}_j$ is calculated by $Q{\text{'}}_j=\frac{\delta {\varphi }_j{}^3}{1-{\varphi }_j}+\left(\frac{1-{\varphi }_j}{{\varphi }_f}\right)\frac{{\varphi }_{j+\frac{1}{2}}-{\varphi }_{j-\frac{1}{2}}}{\frac{1}{2}\Delta {\xi }_{j+1}+\frac{1}{2}\Delta {\xi }_j}$. The average porosity ${\bar{\varphi }}_j$ in the spring element j is given by ${\bar{\varphi }}_j^{k+1}=\frac{V_{cell,j}-V_{cell,dispersed}}{V_{cell,j}}=1-\frac{\left(1-{\varphi }_0\right)\Delta {\xi }_0^{\ast }}{\Delta {\xi }_j^{*k+1}}$. If assuming that the flow rate $Q^{\prime }$ changes approximately linear over each cell, the average porosity represents the middle value of the cell, e.g. ${\varphi }_{j-\frac{1}{2}}={\bar{\varphi }}_j$. We need to interpolate these middle cell porosities at $j-\frac{1}{2}$ to the value at the upstream side cell boundary j, e.g. at the position of “the bead” at ${\xi }_j^{*}$. For cell wall j: ${\varphi }_j^{k+1}={\varphi }_{j-\frac{1}{2}}^{k+1}+\frac{\Delta {\xi }_j^{*k+1}}{\Delta {\xi }_{j+1}^{*k+1}+\Delta {\xi }_{j-1}^{*k+1}}\left({\varphi }_{j+\frac{1}{2}}^{k+1}-{\varphi }_{j-\frac{1}{2}}^{k+1}\right)$. For the lower first cell boundary j = 0: ${\varphi }_0^{k+1}={\varphi }_{\frac{1}{2}}^{k+1}-\frac{\Delta {\xi }_1^{*k+1}}{\Delta {\xi }_1^{*k+1}+\Delta {\xi }_2^{*k+1}}\left({\varphi }_{\frac{3}{2}}^{k+1}-{\varphi }_{\frac{1}{2}}^{k+1}\right)$. Extrapolating for the upper cell upper wall N: ${\varphi }_N={\varphi }_{N-\frac{1}{2}}+\frac{\Delta {\xi }_N^{*}}{\Delta {\xi }_N^{*}+\Delta {\xi }_{N-1}^{*}}\left({\varphi }_{N-\frac{1}{2}}-{\varphi }_{N-\frac{3}{2}}\right)$. 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 ${\displaystyle \underset{j=1}{\overset{N}{\sum }}{\bar{\varphi }}_j\Delta {\xi }_j^{*}}+\left(1-\Theta \right)$, which always equals ${\varphi }_0$, whereby $Q^{\prime }$ is non-negative and Θ increases in each step, otherwise, the simulation stops. An improvement might be to impose the required monotony in $Q^{\prime }$ and $\varphi$, and require the $\frac{\text{d}\varphi }{\text{d}\xi }\le \frac{-\delta {\varphi }^4}{{\left(1-\varphi \right)}^2}$, 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 $U=\frac{\partial \xi }{\partial T}\approx \frac{\partial }{\partial \varphi }\left\{\frac{\delta {\varphi }^3}{1-\varphi }\right\}=\frac{\delta {\varphi }^2\left(3-2\varphi \right)}{{\left(1-\varphi \right)}^2}$. For the cake with $\varphi \le {\varphi }_0$, a maximum rate at $\varphi ={\varphi }_0$ or $U_{\max }=\frac{\delta {\varphi }_0^2\left(3-2{\varphi }_0\right)}{{\left(1-{\varphi }_0\right)}^2}$ occurs. We expect sufficient small displacements to maintain sufficient stability if ΔT = 0.5Δξ/U_max. A Fortran simulation for the bead and spring model for sedimentation WE 6.2% in centrifuge for ${\varphi }_0$ = 0.293 for δ = −0.22, −29, −115, −261 and −464, for N = 10 or 20, appears reproducibly