Particle movements provoke avalanche-like compaction in soft colloid filter cakes

During soft matter filtration, colloids accumulate in a compressible porous cake layer on top of the membrane surface. The void size between the colloids predominantly defines the cake-specific permeation resistance and the corresponding filtration efficiency. While higher fluxes are beneficial for the process efficiency, they compress the cake and increase permeation resistance. However, it is not fully understood how soft particles behave during cake formation and how their compression influences the overall cake properties. This study visualizes the formation and compression process of soft filter cakes in microfluidic model systems. During cake formation, we analyze single-particle movements inside the filter cake voids and how they interact with the whole filter cake morphology. During cake compression, we visualize reversible and irreversible compression and distinguish the two phenomena. Finally, we confirm the compression phenomena by modeling the soft particle filter cake using a CFD-DEM approach. The results underline the importance of considering the compression history when describing the filter cake morphology and its related properties. Thus, this study links single colloid movements and filter cake compression to the overall cake behavior and narrows the gap between single colloid events and the filtration process.


Supplementary files
SI1_ChannelGeometry.stl is the CAD-design file of the microfluidic channel that is used for the filtration experiments.

CFD-DEM Simulation
The study used the simulation framework CFDEM ©. CFDEM © combines computational fluid dynamics (OpenFOAM ©) with a discrete element method (LIGGGHTS ©) to describe particle motion inside a fluid flow. 1, 2 A two way coupling algorithm determines interactions between particle-particle and fluid-particle.

Discrete Element Method
Discrete element method (DEM) is a Lagrangian simulation method developed by Cundall and Strack. 3 The motion and interactions of particles are solved based on Newton's law of motion. Each particle trajectory is solved by the following force balance: 2 where m andẍ i are the mass and the acceleration of the particle, respectively. The motion of each particle depends on contact forces (F n , F t ), body forces (F b ) and forces arising due to the surrounding fluid phase (F f ) . 4 Hertz contact mechanics is used to determine the tangential and the normal contact forces F n and F t : δ n and δt are the normal and tangential overlap of two surfaces. The constants k n and k t are the elastic constants and γ n and γ t the viscoelastic damping constants. The constants k and γ are material properties depending on the Young modulus E, the Poisson ratio ξ and the coefficient of restitution e. 5 The normal and tangential relative velocity is termed v n and v t . The relative tangential velocity includes both the tangential velocity between the spheres and the relative motion due to rotation characterized by the angular velocity.
Besides the Hertzian model, lubrication forces F l ub between the particles are considered in the force balance. If two spheres in a viscous fluid approach each other, the fluid between the spheres has to be squeezed out. The resulting force. When the distance between the surfaces decreases, the pressure gradient necessary to squeeze out the fluid increases. The viscous friction raises through the increased pressure leading to slower relative velocities of the approaching spheres. This lubrication force can be approximated by the following correlation: 6 where η is the fluid viscosity, v rel is the relative velocity between the two surfaces, D is the distance between the two surfaces, a 1 and a 2 are the radii of the two spheres. The lubrication interaction raises to infinity if the separation distance D goes to zero. A cut-off distance D 0 = 10nm is introduced at which the separation distance is assumed to be constant.
A torque balance is applied in addition to the force balance: where I is moment of inertia, ω is the angular velocity and M are the contact torques. This study applied the constant direct torque model to account for rolling friction, 7 which is determined as follows: where C r f defines the coefficient of rolling friction and ω rel,shear is the projection of the relative angular velocity into the shear plane.

Unresolved CFD-DEM approach
The volume averaged CFD-DEM approach is used for modeling the motion of the particle-laden incompressible fluid phase. The volume-averaged Navier-Stokes equations, consisting of the continuity equation and the momentum equation are applied 2 where α f denotes the void fraction of the fluid, u the velocity of the fluid, ρ f the fluid density, p the pressure, R pf the momentum exchange between particle and fluid, and τ = ν(∇u) T the stress tensor with the kinematic viscosity ν. A porous layer is integrated into the CFD part of the simulations representing the pressure drop caused by pore structure in the experiments. The hydrodynamic resistance of the porous layer is added to the Navier-Stokes equation as follows: The resistance k M is chosen to reflect the hydrodynamic resistance of the microfluidic experiments. The momentum exchange R pf is calculated with the following equation: 2 with The momentum exchange is calculated for each CFD cell with the volume V cell , based on the difference between the fluid velocity u and the volume averaged particle velocity v . The drag correlation based on the work of Gidaspow is chosen: 8 The drag pre-factor β p f decreases with the void-fraction ε of the fluid phase: The divided approach is used to calculate the fluid void fraction in a CFD cell. 9 The divided approach splits the particle volume in 29 of equally sized elements. The center point of each element is calculated and the solid fraction is spread over the mesh elements.
Additional smoothing of the void fraction is applied to improve stability of the particle-liquid coupling. 10 The isotropic diffusive smoothing can efficiently be applied to smooth the void fraction field and the momentum exchange term. The following equation is solved for both variables which are represented ξ : The isotropic smoothing method is controlled with the smoothing length λ , which is chosen to be the mean particle's diameter of 20 µm.

Simulation domain and conditions
The simulations were performed in a cuboid including a porous structure, which is shown in Figure 1. The mesh is chosen to be a structured mesh of 25x25x25 µm ensuring the particle size to be smaller than the mesh size. Microgels were filtered on a particle-impermeable layer. The size of the microgels varies over a Gaussian distribution of 18 to 22 micrometers with an average value of 20 micrometer. The particles were filtered with a constant transmembrane pressure (TMP) of 17 mbar until the filter cake reached a height of roughly 250 µm. After filtration, the filter cake is shortly exposed to a TMP of 175 mbar leading to compression of the filter cake. After compression, the cake is relaxed by reducing the TMP again to 17 mbar until steady state is reached. The simulations lasted up to two weeks on eight cores on an AMD Ryzen Threadripper 2990WX processor

Data analysis
The simulation results were analyzed with a python script to gain more information about the cake compression and the cake morphology.
To calculate the filter cake compaction, the filter cake was subdivided into regions of 20 µm in the z-direction. The distance between contacting spheres was calculated under consideration of the periodic boundaries. The resulting contact distances were averaged for each segment, and the standard deviation was calculated. Hence, the degree of compaction can be measured depending on the distance to the membrane.
The coordination number of the particles was calculated to analyze the cake packing. The coordination number describes the number of contacts of each particle. Equal sized spheres in crystalline regions posses in close-packing a coordination number of 12. Again, the particle's location was subdivided into regions of 20 µm in the z-direction. Thereby, the degree of coordination is measured depending on the distance to the membrane.