A note on the modelling of lubrication forces in unresolved simulations

• A rigorously-defined lubrication force model for unresolved simulations is proposed. • Surface roughness and particle deformation are considered. • Numerical stability is ensured by a continuous force formulation. • Good correspondence with experimental data is demonstrated.

• A rigorously-defined lubrication force model for unresolved simulations is proposed.• Surface roughness and particle deformation are considered.• Numerical stability is ensured by a continuous force formulation.• Good correspondence with experimental data is demonstrated.

Keywords:
Lubrication force Fluid-solid systems Multiphase modelling Unresolved simulations CFD-DEM

Introduction
Unresolved simulation methods (e.g.CFD-DEM [1]) are powerful tools to study interacting fluid (gas or liquid) and particle phases.Contrary to fully-resolved simulations (e.g.[2][3][4][5]), in the unresolved methods the fluid flow is only evaluated at scales larger than the particles, relying on models to account for the sub-grid particle-fluid interaction, i.e. the drag force and possibly additional forces [6].It is well-known that the lubrication forces arising during immersed collisions significantly influence the macroscopic behaviour of liquidsolid systems [7,8], and even gas-solid systems with sufficiently small and/or light particles [9].Experimental work on liquid-solid fluidised of the collision as it includes the increased interaction distance between the solid particles, but often relies on arbitrary cut-off distances to ensure numerical stability and/or computational speed.Campos et al. [11] observed increased heterogeneity within the bed, while through a similar approach Wang et al. [12] found a more homogeneous particle distribution resulting from lubrication forces.These contradictory findings illustrate that an accurate description of the lubrication forces and their effects on particle interactions is imperative to simulation of fluid-solid systems.Clearly, a universally-accepted method is currently unavailable, and would greatly benefit the scientific community.Such a method needs to combine accuracy, simplicity, numerical efficiency and robustness to be applicable in large manybody simulations.In this brief communication, we propose a novel implementation of lubrication forces for unresolved simulations, aiming to provide a rigorously defined method to be used in future work.The model is derived based on physical and numerical considerations, and compared with experimental data to demonstrate its accuracy.

Proposed model
The lubrication force is associated with the draining and filling of a thin fluid film between the surfaces of colliding solid bodies.The first solutions of it were presented by Brenner and Cox [13,14] in the 1960s.The normal lubrication force acting on two spherical bodies  and  can be described by Eq. ( 1), where   is the fluid viscosity,   the relative normal velocity,  * =     ∕(  +   ) the reduced radius, ℎ = |⃗   | −   −   the height of the gap between the surfaces, and |⃗   | the vector between the particle centres.
In various experimental studies (e.g.[15]), the rebound of colliding particles in different media was studied.Consistently, a strong dependence of the restitution coefficient  = | ⃗  1 |∕| ⃗  0 | (0 and 1 denoting the approach and rebound, respectively) on the dimensionless Stokes number was observed.Here, we will employ the definition by Yang and Hunt [16] (Eq.(2), where  * =     ∕(  +   ) denotes the reduced mass of the bodies).For collisions below a critical Stokes number St  , no rebound is observed.The exact value of St  is unknown, but values in the range 5 ≤ St  ≤ 20 are commonly found (e.g.[16][17][18]).Conversely, at high Stokes number (St > 2000 [16]), influence of the interstitial medium is negligible.For accurate simulation of fluid-solid systems, it is of vital importance to accurately capture the transition between these two extremes. (2)

Minimum approach
As seen from Eq. ( 1), the lubrication force diverges as ℎ approaches zero.This leads to a paradox where two approaching bodies can never truly come into contact.This paradox relies on the assumption of perfectly rigid and smooth surfaces, which does not hold in practical examples.Joseph and Hunt [19] described that rigid, rough surfaces interact through their asperities, while sufficiently smooth and soft surfaces deform to maintain a finite minimum gap height ℎ min .A solution for the latter case, the elastohydrodynamic collision, was derived by Davis et al. [20].The deformation is governed by the elasticity parameter  (Eq.( 4) [20], where ∕  is the reduced elastic modulus,  is Poisson's ratio, and  ,0 and ℎ 0 denote the initial normal velocity and gap height).A third contact mechanism, glass transition of the liquid phase, was described by Donahue et al. [21].We have not included this mechanism into our model, as it is not significant under the present conditions and introduces the glass transition pressure, a highly challenging parameter to obtain.The approach distance of two rigid surfaces in contact through the tops of their asperities is given by the mean value of their roughness heights  (Eq.( 3), Fig. 1a).The minimum approach due to deformation is expressed by Eq. ( 5) [18] (Fig. 1b).To account for non-monotonous movement of the particles in simulations (e.g.due to interaction with other particles), the highest occurring normal velocity during the approach should be implemented as approach velocity  ,0 .The maximum value of both minimum approach distances is taken as the actual minimum approach distance of two rough, deformable particles (Eq.( 6)).

Normal force
Following the findings of Joseph and Hunt [19], the normal component of the lubrication force can be evaluated independently from the tangential component.The force exerted on two approaching particles is expressed by Eq. ( 8), where the gap height ℎ is truncated to the minimum approach distance ℎ min .This force is applied to any two particles with a separation smaller than a cut-off distance ℎ co , but larger than 0. Here, we used a cut-off distance of ℎ co =  * , as suggested by Yang and Hunt [16].If the two approaching particles overcome the lubrication force and collide, their interaction stemming from the deforming surfaces is described by an established contact model.In this work, the contact force ⃗   is evaluated through the Hertzian contact model, an extensive description of which can be found in the work of Di Renzo and Di Maio [22].In a lubricated collision, the deformation of the surfaces already starts when the minimum approach distance ℎ min is reached.Therefore, the normal overlap   is increased by ℎ min , as shown in Eq. (7).For gap heights between ℎ min and 0, the total force between the particles is a combination of the lubrication force and the contact force, either through asperity contact or elastohydrodynamic deformation.In this transition zone, the total force is approximated by a lever-type expression, which linearly evolves between the lubrication force  , and the contact force  , .This simple expression ensures the total force remains continuous, which is important for application of the model in large-scale simulations, as sudden jumps of force might lead to unstable behaviour during numerical integration, demanding a small integration time step and impeding numerical efficiency.The resulting normal force is summarised by Eq. (9).

Tangential force
Besides the normal force, lubrication effects also play a role in the tangential interaction of bodies.Elaborate analytical expressions are available from O'Neill and Majumdar [23,24], but do not account for the roughness of the surfaces.Instead, we employ empirical approach more suited for large-scale simulations, based on the findings of Joseph and Hunt [19], who extended their experimental work on particle-wall impacts by performing oblique collisions.By observing the rebound angle, they were able to determine the coefficient of friction   .The tangential component of the collision was found to behave much like a dry collision, with a drastically reduced friction coefficient in cases where a thin fluid film remains between the surfaces.It can therefore be modelled through an established contact model describing stickslip behaviour by merely setting the appropriate value for the friction coefficient depending on the contact mechanism (Eq.( 10)).Note that this approach ignores tangential lubrication forces between particle pairs that do not explicitly collide, i.e. which remain separated by more than the minimum approach distance ℎ min .Tangential lubrication forces are relatively weak at such distances, justifying this neglect.The value of the lubricated friction coefficient  lub is not explicitly known.Attempts of Joseph and Hunt [19] to predict the effective coefficient of friction have shown the importance of the local pressure and temperature in the fluid film during the collision.To avoid the need for such elaborate models, the lubricated friction coefficient can be experimentally determined, or estimated.Typically, it can be set one order of magnitude smaller than the dry friction coefficient [19].In this work, we employed the values reported by Yang and Hunt [16].

Performance
For validation, a binary collision model was implemented in MAT-LAB, which is available at [25].The equations of motion were integrated using a velocity Verlet scheme, with a time step equal to 2% of the Hertzian contact time.Collisions corresponding to the work of Yang and Hunt [16] were simulated, and compared with the experimental results.Table 1 lists the used material properties [16]; Fig. 1c defines the geometry of the collisions.Fig. 2 shows the gap height ℎ during collisions of Delrin spheres at various Stokes numbers.For sufficiently low St, there is no rebound and the particles come to a complete halt at a finite gap height.Beyond a critical St  , rebound can be seen to increase with increasing St.
In their experiments of a moving and a stationary sphere, Yang and Hunt [16] reported noticeable movement of the target when ℎ = 0.8 * .Similar behaviour is observed in the numerical results.At St = 1 the spheres never touch, but interact purely hydrodynamically.This is a good indication that the lubrication force increases the range at which particles influence each other, beyond their physical size.Fig. 3 shows the coefficient of restitution  during head-on collisions, comparing the current model with the findings of Yang and Hunt [16].Very similar results were obtained by other groups (e.g.[26,27]).Clearly, the proposed method captures the trends of the experimental data well.Larger deviations are observed at lower Stokers numbers (St < 20), where Yang and Hunt [16] reported strong experimental scatter due to disturbances in the ambient fluid.Above St = 20, the relative error quickly drops to lower values, below 20%.Similar deviations were found for previous models (e.g.[18]), while these models often rely on fitted, estimated or difficultly measured parameters.The current method thus provides satisfactory accuracy, without the need for such precarious values.As expected, smoother materials exhibit a higher critical Stokes number, stemming from their smaller minimum approach distance.The importance of the minimum approach distance is further demonstrated by Fig. 4, which shows the restitution coefficient of Delrin spheres as a function of the surface roughness.The influence of surface roughness is strongest in the transitional regime 20 < St < 1000, where the collision of smoother materials is dampened further.A plateau is reached when ℎ min  = ℎ min  (indicated by the dotted line), when the surfaces deform and the minimum approach distance does not decrease any further.Accounting for this effect greatly improves accuracy, especially for simulation of soft particles and intermediate Stokes numbers.
The model performance during oblique collisions is demonstrated by Fig. 5, which shows impact and rebound angles in three different ranges of the normal Stokes number St  = St ⋅ cos( ,0 ).The lowest St  range (conducted in an 80w% glycerol solution) corresponds to elastohydrodynamic collisions, while the two higher ranges (conducted in water) correspond to asperity contact.Dampening of the collisions leads to decreased rebound angles at lower Stokes numbers.Fig. 5 shows good correspondence between numerical and experimental results, especially for the lowest and highest St  range.Direct comparison for the middle range is precarious, as this range is fairly wide and located around the major transition between viscous and inertial collisions.Considering this, we can only conclude that the current model captures the oblique collisions well.

Closing remarks
In this short communication, we have proposed a lubrication force model for future work on unresolved liquid-solid simulations.We have indicated the importance of the interaction distance and surface deformation for accurate representation of immersed collisions.A comparison with experimental results has shown that the current method captures the physics of these collisions very well, while ensuring numerical stability.Therefore, the proposed method is expected to improve accuracy of existing unresolved fluid-solid models.We encourage the scientific community to review the working example of our model published at [25], and assess its applicability for future work.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Close-up of contact between rigid, rough particles (a), smooth, deformable particles (b) and collision geometry, indicating impact and rebound velocities and angles (c).

Fig. 3 .
Fig. 3. Coefficient of restitution of head-on collisions for various materials, as a function of the Stokes number.Symbols indicate experimental data by Yang and Hunt [16], with error bars indicating data scatter.Lines indicate results from the current model.

Fig. 4 .
Fig. 4. Influence of surface roughness on the restitution coefficient for head-on collisions of Delrin spheres.Dotted line indicates the condition ℎ min  = ℎ min  for increasing Stokes number.

Fig. 5 .
Fig. 5. Angle of rebound  ,1 of oblique collisions for low (St  < 25, Steel I), intermediate (40 < St  < 200, Delrin) and high (500 < St  < 2400, Steel I) normal Stokes numbers, as a function of the impact angle  ,0 .Symbols indicate experimental data by Yang and Hunt [16], lines indicate results from the current model.
Fig. 2. Evolution of gap height ℎ during collisions of Delrin spheres, at various Stokes numbers.